How to Write a Radiative-Convective Model

Barton Paul Levenson
631 Montclair Street
Pittsburgh PA  15217-2808

(412) 521-7040
levenson1960@gmail.com

21,000 words

**1

	Start with a simple, zero-dimensional energy balance model.  Power is energy per unit time.  If the Earth is neither warming nor cooling, then power in must equal power out.  Incoming power is

	PIN = S π R2										|1

where S is the solar constant, π the circle constant (3.14159...) and R the radius of the Earth.  π R2 is Earth's cross-sectional area, treating it as a circle.
	For S = 1360 watts per square meter, and R = 6,371,010 meters, we have PIN = 1.734 x 1017 watts.  Quite a lot of power!

**2

	Earth gets sunlight on its cross-sectional area (π R2), but it radiates from its whole area, which is that of a sphere (4 π R2).  The "flux density" radiated--power per unit surface area--is found from the Stefan-Boltzmann law:

	F = ε σ T4										|2

	Here ε is the "emissivity," or radiative efficiency, which has to be between 0 and 1.  For a planet we assume ε = 1 at the top of the atmosphere.
	σ is the Stefan-Boltzmann constant, 5.670373 x 10-8 W m-2 K-4 in the SI.
	T is the absolute (Kelvin scale) temperature.  So power out is:

	POUT = 4 π R2 σ T4									|3

**3

	In thermal equilibrium--neither heating nor cooling--power in has to equal power out:

	PIN = POUT										|4

	S π R2 = 4 π R2 σ T4								|5

	Terms on both sides cancel, leaving:

	S = 4 σ T4										|6

or

	      S
	T = []0.25									|7
	     4 σ

	This implies a temperature for the Earth, given S = 1360 W m-2, of 278 K.  Earth's actual mean global annual surface temperature is 288 K, so this isn't too far off.  But it's a greatly oversimplified model.

**4

	How can we make this more realistic?  For one thing, we can account for reflection.  Earth isn't solid black; it reflects 30% of the sunlight it gets back out to space.  Some reflection comes from land and sea, more from ice, yet more from clouds.  Some comes from Rayleigh scattering by individual air molecules.  The fraction reflected is the "albedo," A--for climate purposes, the "bolometric Russell-Bond spherical albedo."
	Let's modify our temperature equation appropriately:

	      S (1 - A)
	T = []0.25								|8
	         4 σ

	For S = 1360 W m-2, and A = 0.3, this gives T = 255 K.
	This is below the freezing point of water, 273 K.  So there's a lesson here, first noticed by Jean Joseph-Baptiste Fourier in 1824:  Sunlight alone is not enough to account for the temperature of Earth's surface.  If sunlight were the whole story, Earth should be frozen over!

**5

	Enter the greenhouse effect.
	Fourier knew only that Earth's surface was warmer than sunlight alone accounted for.  He thought the atmosphere must somehow act as a blanket for the "chaleur obscure" (infrared radiation), but he didn't know how.
	John Tyndall began to solve that puzzle with lab work in 1858.  He found that water vapor and carbon dioxide let most sunlight pass through them, but absorbed infrared light ("IR").  They were "greenhouse gases," by a sloppy and very inaccurate analogy with an agricultural greenhouse.
	It works like this.  Sunlight passes through air easily.  It heats the ground.  The warm ground radiates infrared light.
	Greenhouse gases absorb the IR.  They warm up.  They radiate IR of their own--and half of this goes right back down to the ground.  The ground has both sunshine and "airshine" falling on it, so it's warmer than it would be from sunshine alone.

**6

	Does this violate conservation of energy?  Is energy coming out of nowhere in the greenhouse model?
	No, not at all.  Not if you do the math.  So let's do that.
	We will use the next level of modeling from our zero-dimensional model:  a one-dimensional "glass slab model."  It has three parts stacked atop each other:  Space, air, and the ground.
	Sunlight comes out of space in quantity F.  It passes through the air unperturbed.  The ground absorbs it.  The ground radiates IR back up, in quantity G.  The air absorbs all this, and having both a top and a bottom, radiates its own IR, A, both up and down.
	This is where some folks get confused.  Has the air doubled the energy?  No.  Half goes up, half goes down.  The energy balances for each layer follow, with input on the left of the equals sign, output on the right:

1.  Space:		F = A									|9
2.  Air:		G = 2 A								|10
3.  Ground:		F + A = G								|11

	Let's put numbers to it.  Take F as Earth's absorbed climate flux density, (S / 4) times (1 - A).  For solar constant 1360 W m-2 and albedo 0.3, this works out to 238 W m-2.
	Since PIN = POUT (energy must balance for the planet to be in radiative equilibrium), A must equal F (equation 9).  Equation 10 then tells us G is twice this, 476 W m-2.  Assuming perfect emissivity (not far from the truth), the Earth's temperature must then be (476 / σ)0.25 or 302 K.
	We've blown past the observed surface temperature of 288 K.  Our surface has a greenhouse effect of 302 - 255 = 47 K, compared to the actual 288 - 255 = 33 K.  But we're getting the right order of magnitude.  The model doesn't always get more accurate as it gets more realistic (!), but it tends that way.  Generally, the closer you get to an accurate model, the better it will perform.

**7

	Can we improve the model further?
	Sure.  Let's use multiple layers.  Layer 1, at the top, radiates A1 up and down, layer 2 just under it radiates A2 up and down.  In addition to the "G" IR from the ground, air layers can absorb IR from each other.
	It turns out that with two layers, the temperatures are:

	255 K for layer 1.
	302 K for layer 2.
	335 K for the ground.

	For any number of "glass slab" layers, transparent to sunlight but blackbody in the infrared, the temperature of the ground becomes:

	T = Te (N + 1)0.25									|12

where Te is Earth's effective temperature (255 K) and N the number of blackbody layers.  Thus for one layer, T = 302 K, for two layers, T = 335 K, and so on.
	Suppose radiation is the only process involved.  It isn't, but that's a fair approximation.  How many blackbody layers do we need to account for the surface temperatures of Venus, Earth, and Mars?  Solving for N,

	N = (T / Te)4 - 1									|13

	Venus has Te = 184 K and Ts 735 K.*
	Earth has Te = 255 K, Ts = 288 K.
	Mars has Te = 210 K, Ts = 214 K.

*Ts = surface temperature

	So how many blackbody slabs do we need to produce these temperatures?  254 for Venus (!), 0.63 for Earth, 0.078 for Mars.

**8

	How do greenhouse gases absorb IR?  And why only IR?
	It's a quantum effect.  Atoms and molecules move.  They translate (move from place to place), they vibrate, they rotate.  But due to quantum effects, they can only vibrate and rotate in multiples of basic values.  Because of this, radiatively active molecules only absorb or emit at particular "absorption" lines in the spectrum.
	These lines are all through the spectrum.  Each atom and molecule has its particular "fingerprint" of lines.  That's how we know what the sun and stars are made of, a feat they thought, in the early 19th century, would always be impossible.  For greenhouse gases, there are few absorption lines in the visual, but many in the infrared.  That's not too much of a surprise; the infrared range is a lot bigger than the visual range.

**9

	A perfect or "blackbody" radiator distributes its radiation by the rather complicated Planck function:

	               c1
	I(λ) = ——————————————————							|14
	       λ5 exp(c2 / [λ T])

where c1 is the "first radiation constant," c2 the "second radiation constant," λ the wavelength, T temperature, and I "irradiance" in watts per square meter, per meter of wavelength, per steradian of solid angle.
	If you graph the function, it looks like a bell curve leaning over a bit to the left.  It drops to zero at wavelength 0, rises to a peak, then gradually drops toward zero again as wavelength heads out to infinity.  But though the curve is infinitely long, the area under it is finite.  If you've taken calculus you've run into this sort of thing before.
	The total amount radiated, the integral of the function over all wavelengths, is

	P = σ T4										|15

which is the Stefan-Boltzmann law.

**10

	The location of the Planck function's peak depends only on temperature.  You find it by the very simple "Wien's Law:"

	      c3
	λmax = ——										|16
	       T

where c3 is the "third radiation constant" and T the absolute temperature.  Wavelengths of light are most conveniently measured in microns (1 μ = 10-6 meter), for which c3 = 2897.7721 μ K.  Thus for the sun, with Te = 5774 K, peak radiation is at 0.5 μ.
	Actually, the Sun isn't a perfect blackbody, and the peak actually comes at about 0.48 μ.  Sorry about that.
	By way of contrast, the Earth's surface, at 288 K, has peak radiation at 10 microns, in the IR range.  That's why the sun gives off visible light and the Earth's surface gives off infrared light, and that's why the greenhouse effect works.

**11

	The most accurate radiative transfer models use "line-by-line" calculations, where the absorption from each line of all the greenhouse agents involved is calculated, crawling across the spectrum (mathematically) to handle 100,000 or so absorption lines.  These models are agonizingly slow, so quicker approximations are always welcome.
	One way is to chart where the lines cluster thickly together and consider them "bands" with a certain mean absorption.  There are many kinds of band model, involving much mathematics and hundreds, if not thousands, of published papers.  It's simpler still to just grab someone else's band model from an earlier paper, which is what I usually do when writing my "radiative-convective" models.

**12

	A simple model of light transmission through a medium is "Beer's Law" (actually the Beer-Lambert-Bouguer law, but who wants to remember that?).  Given certain circumstances (monochromatic light, no scattering, homogeneous medium, etc.), the amount of light that gets though a slab of medium is

	T = exp(-k ρ ds)									|17

where T is the transmittance (FOUT over FIN), a number between 0 and 1.  k is an "absorption coefficient," ρ the density of the absorber, and ds the distance traveled.
	Technically, k applies only to a given wavelength (or frequency), but in a band model, you take an average over a whole range of wavelengths (or frequencies).

**13

	Ignoring chemical reactions or phase changes, light can only have three fates when it hits a medium:  it can be absorbed, reflected, or transmitted.  In a three-dimensional setup, we say "scattered" instead of reflected, but you get the idea.  The fractions absorbed, reflected, and transmitted then have to sum to 100%:

	A + R + T = 1									|18

	This is "the radiation identity."

**14

	Let's do a radiative transfer problem.  How much infrared light between 14 and 16 microns wavelength gets through a slab of air one meter deep?  Let's assume carbon dioxide is the only absorber present, at 400 parts per million by volume.
	The CO2 absorption coefficient in this range is 20 square meters per kilogram--it's an "absorption cross-section," like the cross-sections nuclear physicists use for atomic reactions.  The optical depth, tau, is:

	τ = k ρ ds										|19

	k here is 20 m2 kg-1, ds is 1 m.  ρ is harder.  Let's assume the air is at room temperature, 293 K, and sea-level pressure.  These conditions are called "Standard Temperature and Pressure," or STP.
	The ideal gas law is:

	P V = n R T										|20

where P is pressure, V volume, n the number of particles (atoms, molecules, whatever), R the universal gas constant, and T absolute temperature.  For STP, P = 101,325 pascals and T = 293 K.  Our volume will be 1 m3.  R, the universal gas constant, has the value 8314.4621 joules per kilomole per kelvin (J kmol-1 K-1).
	A kilomole is how many particles are present when you've got the atomic weight (or molecular weight) of the substance in kilograms:  About 6.02214 x 1026.  Solving for n, we find that a cubic meter of air at STP has this many molecules in it:

	                  101,325 x 1
	n = P V / R T = ——————————————— = 0.0416 kmol				|21
	                8314.4621 x 293

	The volume fraction is also the molar fraction, so 0.0004 of this--0.0000166 kmol--is carbon dioxide.  The molecular weight of CO2 is 44.0098 amu, so we've got 0.00073 kilograms of CO2 in our slab--not a lot!
	The optical depth, then, is τ = 20 x 0.00073 x 1 = 0.0146.  The transmittance is the exponential of the negative of this, T = 
exp(-0.0146) = 0.986.  98.6% of the light got through the slab, 1.4% was absorbed by CO2.  So far so good.

**15

	I introduced optical depth rather hastily in that last post.  Again, the mathematical definition:

	τ = k ρ ds										|22

where k is absorption coefficient, ρ density of the absorber, and ds the path length.
	Note the units.  If k is in m2 kg-1, ρ in kg m-3, and ds in m, then τ comes out dimensionless.  But what is "optical depth?"
	It's a measure of resistance to the passage of light.  A high-optical-depth medium lets very little light through, a low-optical-depth medium is more transparent.  The definition is such that an optical depth of τ = 1 diminishes the passage of light by a factor of Euler's number, e = 2.71828..., the base of the natural logarithms.
	Media with τ << 1 are "optically thin."  Those with τ >> 1 are "optically thick."  Thus Venus has an optically thick atmosphere, while Mars has an optically thin one (except during dust storms!).

**16

	The absorption coefficients, pressure, temperature, and even composition of air change with altitude.  Thus came about the idea of the "column model" of a climate system.  You divide the atmosphere into layers, and add a layer of ground.  The air layers can all be the same thickness, or different thicknesses, depending on what you want to find out.  A model like this treating radiation only is a "radiative equilibrium model" or REM.  A more sophisticated treatment, taking into account cooling of the surface at the expense of the atmosphere, is a "radiative-convective model" or RCM.  REMs go way back.  RCMs were first proposed by Hulbert in 1931, and first successfully done on a computer by Manabe and Strickler in 1964.
	A modern RCM, whether of Earth or of some other world (Venus, Mars, Titan, an exoplanet...) usually has dozens or even hundreds of layers.  I'll build a model of Earth's mean climate here, using twenty layers of air and one of ground.  The language will be Visual Basic 2010 Express.
	I'll keep it simple:  console I/O only, no classes or objects.  If you can follow along, this shouldn't be hard to duplicate in your language of choice.  You will need to know some computer programming, though.

**17

	We'll write an REM and gradually turn it into an RCM.  Not professional quality, but pointing the way to how to write a pro-quality model.
	The most important data structures we'll need are absorbers, bands, and layers.  Layers are easiest to visualize, so I'll start with those.
	There are five meteorological layers to the atmosphere.  From the ground up, these are the troposphere, stratosphere, mesosphere, thermosphere, and exosphere.  For our model, we can ignore the last three, since almost all the atmosphere's mass is in the lowest two levels--80% troposphere and 20% stratosphere.
	The lower (troposphere) layer is where we live.  The name is from Greek for "turning," because it's the level where weather takes place--winds, inversions, hurricanes, etc., etc.
	The stratosphere, above that, is from Greek for "layered," because they once had the idea that there was no weather there, no turbulence, and the gases lay atop one another ordered by density.  The stratosphere really does have very little convection, and is much calmer than the troposphere, but the gases in it are still well-mixed.  You don't get layering by molecular weight until very, very high in the very, very thin atmosphere.

**18

	Trenberth et al. (2005) estimated the mass of Earth's atmosphere as 5.148 x 1018 kg.  It would be closer to 5.270 x 1018 kg, estimating from sea-level pressure, but continents and islands stick up above sea-level and take up some of the volume.  Of the atmosphere mass, an average of 1.27 x 1016 kg is water vapor, the most important greenhouse gas.
	We'll divide the air layers evenly by mass, considering the top four as the stratosphere and the sixteen below that the troposphere.  For the ground, we'll use one meter of seawater--a "swamp" model.  This is a common tactic with RCMs.  For the even more elaborate GCMs (General Circulation Models, or Global Climate Models), they usually use many levels of ocean, not to mention rock, soil, ice, cities, etc., but that's because they're trying to track what happens to the whole climate system over time.  RCMs are more for getting a snapshot of the situation at one moment.
	The average density of seawater is 1,025 kg m-3, and its "specific heat capacity at constant pressure" is 3,890 J K-1 kg-1.  That means you need 3890 joules to raise the temperature of one kilogram of the stuff by one kelvin (or degree Celsius, which is identical to a kelvin).
	The sea-level density of air is 1.225 kg m-3 at mean global conditions.  Unlike water, air is compressible, so the density varies markedly with altitude.  We'll ignore density and use pressure to mark layers.  This means we need a different definition of optical thickness:

	τ = k SM										|23

where k is the absorption coefficient and SM is "specific mass" or "mass path"--the mass per unit area.  Note that this has the same units as density times length in the traditional definition of τ.
	The specific heat of dry air is about 1,004 J K-1 kg-1.  The figure for water vapor is higher, about 1,930.  But since water vapor is such a minor constituent, usually about 0.4% by volume, it doesn't much affect the air specific heat.  We can get fairly good answers just assuming a fixed value of 1,006.  And mean molecular weight, which would be 28.9644 for dry air, becomes about 28.94 for normally wet air, since water vapor is fairly light at 18.0153 amu per molecule.

**19

	There are many absorbers in Earth's atmosphere, the most important being water vapor, carbon dioxide, ozone, methane, and nitrous oxide.  But for a good first approximation, we can use water vapor and CO2 alone.  CO2, like CH4 and N2O, is a "well-mixed gas" through 99% of the atmosphere.  But water vapor and ozone vary all over the place, with local conditions as well as with altitude.
	You can get a pretty good global mean for water vapor using its "relative humidity:"

	RH = pH2O / eS									|24

where pH2O is actual water vapor pressure, and eS is "saturation vapor pressure."  If you assume fixed relative humidity at a given pressure level, you get a fairly accurate vertical distribution for water vapor.  It turns out to be very shallow compared to the mixed gases.  Since water is a "volatile" at Earth conditions--it changes phases easily--it stays near the oceans, trapped in that cycle we all learned about in fifth grade--evaporation, condensation, precipitation.  An average molecule of water vapor stays in the air nine days, a very fast cycle indeed.
	Manabe and Strickler used fixed "absolute humidity" in their 1964 model--a given mass of water vapor in a given layer.  As it turned out, it was more realistic to use fixed relative humidity, and Manabe and Wetherald did just that in an improved model in 1967.  Their algorithm was:

	         p - 0.02
	RH = RH0 —————————								|25
	         1 - 0.02

where RH0 is sea-level relative humidity and p is pressure in atmospheres.  Of course, this gives negative humidity when the pressure drops below 0.02 atm, which is physically impossible.  In that case, they just substituted the minimum observed level of RH, which occurs in the stratosphere at about RH = 0.000005.  The figure they used for RH0 was 77%, and most modern RCMs still use that, although you can find figures from 75 to 80%.

**20

	Given all that, here's the first look at our new model:  A table of properties of the layers.
 
		————————————————————————————————————————

		L#	 m (kg)	 cp	p (Pa)	 RH	

		————————————————————————————————————————

		 1	2.574E+17	1006	  2474	0.004	
		 2	2.574E+17	1006	  7423	0.043	
		 3	2.574E+17	1006	 12372	0.083	
		 4	2.574E+17	1006	 17321	0.123	
		 5	2.574E+17	1006	 22270	0.163	
		 6	2.574E+17	1006	 27219	0.203	
		 7	2.574E+17	1006	 32167	0.243	
		 8	2.574E+17	1006	 37116	0.283	
		 9	2.574E+17	1006	 42065	0.323	
		10	2.574E+17	1006	 47014	0.362	
		11	2.574E+17	1006	 51963	0.402	
		12	2.574E+17	1006	 56912	0.442	
		13	2.574E+17	1006	 61860	0.482	
		14	2.574E+17	1006	 66809	0.522	
		15	2.574E+17	1006	 71758	0.562	
		16	2.574E+17	1006	 76707	0.602	
		17	2.574E+17	1006	 81656	0.642	
		18	2.574E+17	1006	 86605	0.681	
		19	2.574E+17	1006	 91553	0.721	
		20	2.574E+17	1006	 96502	0.761	
		21	5.228E+17	3890	106351	1.000
		————————————————————————————————————————

	To calculate the mean pressure in our swamp ocean, I used the engineering standard g:  sea-level gravity of 9.80665 m s-2, and an Earth surface area of 5.10066 x 1014 square meters.

**21

	Let's get coding.  Create a console app named rE in Visual Studio--or the equivalent in whatever you've got.  Rename the main module rEmain.  We'll stick to one module and one source file to keep things simple.
	Import these libraries at the top:

Imports System.Console
Imports System.Math

	Create a new type above the main module.  The layer type may look incredibly complicated.  Don't panic.  We won't use everything in it right away.  But I thought it was better to lay the whole thing out at once, so we have all the parts when we need them later.

    Public Structure layerType
        Dim A() As Double, R() As Double, T() As Double
        Dim Acld() As Double, Rcld() As Double, Tcld() As Double

        Dim Av As Double, Rv As Double, Tv As Double
        Dim Ai As Double, Ri As Double, Ti As Double

        Dim AvCld As Double, RvCld As Double, TvCld As Double
        Dim AiCld As Double, RiCld As Double, TiCld As Double

        Dim above() As Double
        Dim below() As Double
        Dim radOut() As Double
        Dim refUp() As Double
        Dim refDown() As Double

        Dim cp As Double
        Dim cpm As Double
        Dim fcld As Double
        Dim fclr As Double
        Dim gamma As Double
        Dim height As Double
        Dim mass As Double
        Dim mfrac As Double
        Dim MW As Double
        Dim num As Double
        Dim oldTemp As Double
        Dim p As Double
        Dim pAtm As Double
        Dim RH As Double
        Dim temp As Double
        Dim z As Double

        Dim SM() As Double
    End Structure ' layerType

	I'll go over the parts a bit at a time.  Again, don't panic--we won't be using most of this for a long time.
	A, R, and T are the layer's absorptance, reflectance, and transmittance in each band.  These values apply to clear air.  They will be the only radiative properties we use at first.
	Acld, Rcld, and Tcld are the same type of array, but for clouds.
	Av, Rv, and Tv are the mean values of A, R, and T in the "visual" (shortwave) domain.  Ai, Ri, and Ti are the same for the infrared (longwave). 
	AvCld through TiCld are the same values for clouds.
	The arrays above and below refer to the power entering a layer from above and below, respectively.  radOut is the power emitted, and will usually be the same in both directions.  The exception is the ground, which doesn't radiate downward.  And refUp and refDown are the power reflected away in the up and down directions.
	cp is the layer's specific heat capacity at constant pressure.
	cpm is the heat capacity times the mass.  It's there so we don't have to repeat this calculation over and over as the program runs.
	fcld is the fraction of each layer covered by clouds.  fclr is the clear fraction.  The two are related by

	fclr + fcld = 1									|26

	gamma is the layer's lapse rate (heat decrease per unit height, measured in K km-1).
	height is the layer's height in meters.
	mass is its mass in kilograms.  mfrac is the mass fraction compared to the whole atmosphere.
	MW is the layer's mean molecular weight in u (atomic mass units).
	num is the number of kilomoles of molecules in the layer.  mass, MW, and num are also related:

	num = mass / MW									|27

	oldTemp is the layer's temperature on the previous cycle of the program.
	p is the layer's mean pressure in pascals.  pAtm is the pressure in atmospheres:

	pATM = p / P0									|28

here P0 is the standard sea-level pressure, 101,325 Pa.
	RH is the layer's relative humidity, as a decimal.
	temp is its temperature (K).
	z is its mean altitude in meters.
	Finally, SM is the specific mass, or mass path, for each constituent.

**22

	We will need a lot of global constants for various equations, and to size arrays properly.  I'll introduce them a section at a time.

    ' absorber
    Public Const Nabs As Integer = 10
    Public Const Ngases As Integer = 8
    Public Const Ar As Integer = 1
    Public Const CH4 As Integer = 2
    Public Const CO2 As Integer = 3
    Public Const H2O As Integer = 4
    Public Const N2 As Integer = 5
    Public Const N2O As Integer = 6
    Public Const O2 As Integer = 7
    Public Const O3 As Integer = 8
    Public Const ice As Integer = 9
    Public Const wtr As Integer = 10

	Here Nabs is the total number of absorbers (or potential absorbers) in the atmosphere.  Ngases is the total number of gases--the absorbers that aren't gases are cloud particles.  The last ten constants are "symbolic constants"--handy placeholders, because it's easier to understand SM(CO2) as the specific mass of carbon dioxide than it is to remember that the 3 in SM(3) means carbon dioxide.

    ' band
    Public Const Nbands As Integer = 45
    Public Const swlo As Integer = 1
    Public Const swhi As Integer = 19
    Public Const lwlo As Integer = swhi + 1
    Public Const lwhi As Integer = Nbands
    Public Const Nsw As Integer = swhi - swlo + 1
    Public Const Nlw As Integer = lwhi - lwlo + 1

	Nbands, obviously, is the number of bands in our band scheme.  swlo and swhi are the boundaries of the shortwave region, lwlo and lwhi the boundaries of the longwave region.  Finally, Nsw and Nlw are the number of bands in each region.

    ' cloud
    Public Const cldHi As Integer = 6
    Public Const cldMid As Integer = 12
    Public Const cldLo As Integer = 17
    Public Const NChou As Integer = 5
    Public Const ri As Double = 30.0#
    Public Const rw As Double = 5.0#

	cldHi, cldMid, and cldLo are the levels our high, middle, and low cloud decks will be located.  NChou is the number of Chou et al. cloud scheme bands, which will be discussed later.  ri and rw are the mean effective radii of ice and water cloud particles, respectively, in microns.

    ' Earth
    Public Const area As Double = 510066000000000.0#
    Public Const g As Double = 9.80665#
    Public Const mAtm As Double = 5.148E+18#
    Public Const P0 As Double = 101325.0#
    Public Const S As Double = 1361.0#
    Public Const T0 As Double = 288.15#
    Public Const TOA As Double = 0.25# * S

	These are properties of the Earth, obviously--and these are the constants that will change the most if we do an RCM of another plant.  area is the surface area in square meters.  g is mean surface gravity in m s-2.  mAtm is the atmosphere mass in kilograms.  P0 and T0 are the reference (sea-level) pressure and temperature.  S is the Solar constant in W m-2, and TOA is "top of atmosphere" mean illumination, or S / 4.

    ' layer
    Public Const N As Integer = 20
    Public Const all As Integer = N + 1
    Public Const ground As Integer = all

	We have N = 20 layers of atmosphere, and N + 1 or 21 layers including the layer of ground at the bottom.  Finally, ground is another symbolic constant to make statements in code more readable.

    ' physical
    Public Const Rstar As Double = 8314.4621#
    Public Const sb As Double = 0.00000005670373#

	These are physical constants from the sciences.  Rstar (R*) is the universal gas constant, in J kmol-1 K-1.  sb (σ) is the Stefan-Boltzmann constant, in W m-2 K-4.
	That's it for the constants.

**23

	Now, add one global variable, a big array:

    Public layer(All) As layerType

	Now we can write active code.  We'll start with set-up routines:

    ' setPhysAir sets air layer physical properties.
    Sub setPhysAir()
        Const layerMass As Double = mAtm / N
        Const RH0 As Double = 0.8#
        Const kRH As Double = RH0 / 0.98#
        Const minRH As Double = 0.000005#

        For L As Integer = 1 To N
            With layer(L)
                .cp = 1006.0#
                .gamma = 0.0#
                .height = 0.0#
                .mass = layerMass
                .mfrac = .mass / mAtm
                .MW = 28.94#
                .num = .mass / .MW
                .oldTemp = 137.0#
                .p = (L - 0.5#) * layerMass * g / area
                .pAtm = .p / P0
                .temp = 255.0#
                .z = 0.0#

                .cpm = .cp * .mass
                .RH = kRH * (.pAtm - 0.02#)
                If .RH < minRH Then
                    .RH = minRH
                End If
            End With
        Next
    End Sub ' setPhysAir

    ' setPhysGround sets ground physical properties.
    Sub setPhysGround()
        With layer(ground)
            .cp = 3890.0#
            .gamma = 0.0#
            .height = 1.0#
            .mass = 1025.0# * area
            .mfrac = .mass / mAtm
            .MW = 18.02#
            .num = .mass / .MW
            .oldTemp = 137.0#
            .p = 106351.0#
            .pAtm = .p / P0
            .RH = 1.0#
            .temp = 255.0#
            .z = -0.5#

            .cpm = .cp * .mass
        End With
    End Sub ' setPhysGround

    ' setup runs all the set-up routines.
    Sub setup()
        setPhysAir()
        setPhysGround()
    End Sub

	Call setup from the main routine.  We now have a working, although not very useful, program.  Run it.  If there are any errors, debug it.

**24

	The program can be made more interesting by adding an I/O utility and a display routine.  The display routine may be trivial, but it will be just one of many diagnostic screens that can help us find problems:

    ' pause lets the user read before going on.
    Sub pause()
        SetCursorPosition(0, 23)
        Write("Press [Enter] to go on.  ")
        Dim st As String = ReadLine()
    End Sub ' pause

    ' showPhys shows layer physical properties.
    Sub showPhys()
        Clear()
        WriteLine("L#    m (kg)     cp   P (Pa)   RH    T (K)")
        WriteLine("--   ---------  ----  ------  -----  -----")

        For L As Integer = 1 To ground
            Write(L.ToString.PadLeft(2))

            With layer(L)
                Write(.mass.ToString("e2").PadLeft(12))
                Write(.cp.ToString("f0").PadLeft(6))
                Write(.p.ToString("f0").PadLeft(8))
                Write(.RH.ToString("f3").PadLeft(7))
                Write(.temp.ToString("f1").PadLeft(7))
            End With

            WriteLine()
        Next

        pause()
    End Sub ' showPhys

	If you're not doing so already, I recommend putting subroutines and functions in alphabetical order by name.  It makes them much easier to find later while debugging.

**25

	Call showPhys right after the call to setup in the main routine.  Then you should see a console screen appear:

L#    m (kg)     cp   P (Pa)   RH    T (K)
--   ---------  ----  ------  -----  -----
 1   2.57e+017  1006    2474  0.004  255.0
 2   2.57e+017  1006    7423  0.043  255.0
 3   2.57e+017  1006   12372  0.083  255.0
 4   2.57e+017  1006   17321  0.123  255.0
 5   2.57e+017  1006   22270  0.163  255.0
 6   2.57e+017  1006   27219  0.203  255.0
 7   2.57e+017  1006   32167  0.243  255.0
 8   2.57e+017  1006   37116  0.283  255.0
 9   2.57e+017  1006   42065  0.323  255.0
10   2.57e+017  1006   47014  0.362  255.0
11   2.57e+017  1006   51963  0.402  255.0
12   2.57e+017  1006   56912  0.442  255.0
13   2.57e+017  1006   61860  0.482  255.0
14   2.57e+017  1006   66809  0.522  255.0
15   2.57e+017  1006   71758  0.562  255.0
16   2.57e+017  1006   76707  0.602  255.0
17   2.57e+017  1006   81656  0.642  255.0
18   2.57e+017  1006   86605  0.681  255.0
19   2.57e+017  1006   91553  0.721  255.0
20   2.57e+017  1006   96502  0.761  255.0
21   5.23e+017  3890  106351  1.000  255.0
Press [Enter] to go on.

**26

	Let's add another structure:

    Public Structure absorberType
        Dim cp As Double    ' Specific heat
        Dim MW As Double    ' Molecular weight
        Dim name As String  ' Name (formula)
        Dim q As Double     ' Mass fraction.
    End Structure

	For the moment we'll restrict "absorbers" to gases.  Later we can add cloud particles.

**27

	You can get a pretty close match to standard air with just eight gases:  the big ones, nitrogen, oxygen, and argon, none of which are very radiatively active, and five trace gases:  water vapor, carbon dioxide, ozone, methane, nitrous oxide.  The latter five are all greenhouse gases, but we won't use the last three for now.
	Most breakdowns of gases in the atmosphere give the volume fraction.  To get the mass fraction, which is more convenient for RCMs, you need to multiply by the ratio of the individual gas's molecular weight to the mean weight of the medium it's in:

	        μi
	qi = Xi ——										29
	        μ0

where q is the mass fraction for gas i, X is the volume fraction, and μ the molecular weight in AMU.  That gives us:

gas	
cp	
MW	
X	
q

Ar	
 520	
39.948	
0.00934	
0.0129	
CH4	2220	16.0428	0.0000016	0.00000089	
CO2	 844	44.0098	0.00032	0.00049	
H2O	1930	18.0153	0	0	
N2	1040	28.0134	0.78084	0.75522	
N2O	 880	44.0129	0.00000032	0.00000049	
O2	 919	31.9988	0.20946	0.23141	
O3	 818	47.9982	0	0

	This results in a mean molecular weight of 28.9637, darned close to the canonical value of 28.9644.  Note that I'm assuming zero fractions for water vapor and ozone; they'll be added separately.

**28

	Let's add another big array:

    Public absorber(Nabs) As absorberType

	And some setup routines:

    ' set1absorber sets an absorber's properties.
    Sub set1absorber(ByVal a As Integer, ByVal cc As Double, ByVal mm As Double, _
    ByVal nn As String, ByVal qq As Double)
        With absorber(a)
            .cp = cc
            .MW = mm
            .name = nn
            .q = qq
        End With
    End Sub ' set1absorber

    ' setAbsorbers sets gases and cloud particles.
    Sub setAbsorbers()
        set1absorber(Ar, 520.0#, 39.948#, "Ar", 0.0129#)
        set1absorber(CH4, 2220.0#, 16.0428#, "CH4", 0.00000089#)
        set1absorber(CO2, 844.0#, 44.0098#, "CO2", 0.00049#)
        set1absorber(H2O, 1930.0#, 18.0153#, "H2O", 0.0#)

        set1absorber(N2, 1040.0#, 28.0134#, "N2", 0.75522#)
        set1absorber(N2O, 880.0#, 44.0129#, "N2O", 0.00000049#)
        set1absorber(O2, 919.0#, 31.9988#, "O2", 0.23141#)
        set1absorber(O3, 818.0#, 47.9982#, "O3", 0.0#)
    End Sub ' setAbsorbers

**29

	The last (and simplest) new data structure we need is the band type:

    Public Structure bandType
        Dim lo As Double        ' Low boundary.
        Dim hi As Double        ' High boundary.
        Dim solFrac As Double   ' Solar radiation fraction.
    End Structure ' bandType

	We will need to modify some of these data structures later, mostly by adding new components.   But for now, these will do.  Let's add a new big array, and two more setup routines:

    Public band(Nbands) As bandType

    ' set1band sets properties for a single band.
    Sub set1band(ByVal b As Integer, ByVal ll As Double, ByVal hh As Double)
        With band(b)
            .lo = ll
            .hi = hh
            .solFrac = 0.0#
        End With
    End Sub ' set1band

    ' setBands creates the band scheme.
    Sub setBands()
        set1band(1, 0.175#, 0.225#)
        set1band(2, 0.225#, 0.245#)
        set1band(3, 0.245#, 0.28#)
        set1band(4, 0.28#, 0.295#)
        set1band(5, 0.295#, 0.31#)

        set1band(6, 0.31#, 0.32#)
        set1band(7, 0.32#, 0.4#)
        set1band(8, 0.4#, 0.7#)
        set1band(9, 0.7#, 1.8#)
        set1band(10, 1.8#, 1.9#)

        set1band(11, 1.9#, 2.0#)
        set1band(12, 2.0#, 2.1#)
        set1band(13, 2.1#, 2.2#)
        set1band(14, 2.2#, 2.5#)
        set1band(15, 2.5#, 2.6#)

        set1band(16, 2.6#, 2.9#)
        set1band(17, 2.9#, 3.0#)
        set1band(18, 3.0#, 3.5#)
        set1band(19, 3.5#, 4.1#)
        set1band(20, 4.1#, 4.5#)

        set1band(21, 4.5#, 5.0#)
        set1band(22, 5.0#, 6.5#)
        set1band(23, 6.5#, 7.0#)
        set1band(24, 7.0#, 8.0#)
        set1band(25, 8.0#, 9.075#)

        set1band(26, 9.075#, 9.265#)
        set1band(27, 9.265#, 9.32#)
        set1band(28, 9.32#, 9.4#)
        set1band(29, 9.4#, 9.475#)
        set1band(30, 9.475#, 9.61#)

        set1band(31, 9.61#, 9.745#)
        set1band(32, 9.745#, 9.855#)
        set1band(33, 9.855#, 9.95#)
        set1band(34, 9.95#, 10.045#)
        set1band(35, 10.045#, 10.135#)

        set1band(36, 10.135#, 12.0#)
        set1band(37, 12.0#, 13.235#)
        set1band(38, 13.235#, 13.545#)
        set1band(39, 13.545#, 14.0#)
        set1band(40, 14.0#, 16.0#)

        set1band(41, 16.0#, 16.66#)
        set1band(42, 16.66#, 25.0#)
        set1band(43, 25.0#, 30.0#)
        set1band(44, 30.0#, 40.0#)
        set1band(45, 40.0#, 125.0#)
    End Sub ' setBands

	Looking at the Planck curves for Solar and terrestrial radiation, 99% of solar radiation falls in the range below 4 microns, while 99% of the terrestrial falls above it.  Thus the first 19 of our 45 bands are "shortwave," the rest "longwave."
	Why these band boundaries?  They correspond to absorption bands we'll use later for our greenhouse gases.  Until then, just consider them arbitrary.

**30

	The Solar fraction in each band is not, of course, zero, as I have it so far.  But finding the fraction of solar radiation between two wavelengths is one tricky business.  The Planck curve, as noted earlier, is complicated.  Worse, it falls into a class of equations that can't be integrated.  The best you can do is approximate it.
	Fortunately, they worked out such approximations early in the 20th century.  These are infinite-series approximations, but since each extra term is smaller than the last, the series converge.  20 terms generally does it, which is easy on a computer.  I use 200 just to be sure.

    ' D returns the Planck fraction from 0 to a given wavelength.
    Function D(ByVal lambda As Double, ByVal T As Double) As Double
        Const c2 As Double = 14387.77#
        Const k As Double = 15.0# / (PI ^ 4)

        Dim q As Double = c2 / (lambda * T)
        Dim sum As Double = 0.0#

        For m As Integer = 1 To 200
            Dim f As Double = m * q
            Dim term1 As Double = Exp(-f) / (m ^ 4)
            Dim term2 As Double = f * (f * (f + 3.0#) + 6.0#) + 6.0#
            sum += (term1 * term2)
        Next m

        Return k * sum
    End Function ' D

**31

	So modify that one line in set1band to read:

	.solFrac = D(hh, 5774.0#) - D(ll, 5774.0#)

	I'm using 5,774 K for the effective temperature of the Sun.

**32

	Why am I throwing more and more code into this apparently useless program?
	Because when I have enough of the parts working, I can use it to run a "time-marching" algorithm that will predict the temperature of Earth's surface and atmosphere.  Here's how it works:

1.	Start off with an "isothermal" system--all layers the same temperature.  I use Earth's effective temperature, 255 K.
2.	Calculate how much radiation enters and leaves each layer, from the sun and from other layers.
3.	Find the temperature change in each layer as a result.
4.	If no layer has changed temperature by more than a tiny bit, quit.
5.	Otherwise, go back to step 2.

	For an RCM, as opposed to an REM, there's more to it--you have to check after #3 that no layer exceeds a critical lapse rate, and simulate convection if it has.  But we'll get to that later.

**33

	Sadly, VB won't let you declare fixed-length arrays in a structure, which is why we declared all the arrays in layerType with empty parentheses.  Before we get to radiative properties, therefore, we will have to set array lengths.  Lots of things blow up if you try using undimensioned arrays, so this routine should be the first called from setup:

    ' dimension sets array sizes.
    Sub dimension()
        For L As Integer = 1 To ground
            With layer(L)
                ReDim .A(Nbands)
                ReDim .R(Nbands)
                ReDim .T(Nbands)

                ReDim .Acld(Nbands)
                ReDim .Rcld(Nbands)
                ReDim .Tcld(Nbands)

                ReDim .above(Nbands)
                ReDim .below(Nbands)
                ReDim .radOut(Nbands)
                ReDim .refUp(Nbands)
                ReDim .refDown(Nbands)

                ReDim .SM(Nabs)
            End With
        Next
    End Sub ' dimension

**34

	Eventually, we want to set A, R, T etc. by real physical processes.  But it's complicated, so we'll start by simply assigning plausible values.  Here are the setup routines:

    ' setRadAir sets air layer radiative properties.
    Sub setRadAir()
        For L As Integer = 1 To N
            With layer(L)
                For b As Integer = swlo To swhi
                    .A(b) = 0.0#
                    .R(b) = 0.003#
                    .T(b) = 0.997#
                Next b

                For b As Integer = lwlo To lwhi
                    .A(b) = 0.15#
                    .R(b) = 0.0#
                    .T(b) = 0.85#
                Next b
            End With
        Next L
    End Sub ' setRadAir

    ' setRadGround sets ground radiative properties.
    Sub setRadGround()
        With layer(ground)
            For b As Integer = swlo To swhi
                .A(b) = 0.86#
                .R(b) = 0.14#
                .T(b) = 0.0#
            Next b

            For b As Integer = lwlo To lwhi
                .A(b) = 0.996#
                .R(b) = 0.004#
                .T(b) = 0.0#
            Next b
        End With
    End Sub ' setRadGround

	There is almost no IR scattering in Earth's atmosphere.  Visual scattering for air layers is set by taking the observation that Rayleigh scattering reflects about 6% of the incoming sunlight; I just divided it among the air layers by mass.  For the ground, transmissivity must be 0 everywhere (the ground is opaque).  The visual reflectivity is from Project CERES observations; the IR reflectivity is an estimate from the ECHAM5 GCM.  Values for air IR absorptance is after a suggestion by Prof.  Lawrence Dempsey of UC Santa Barbara, who set it at 15% for twenty equal-mass layers.

**35

	With 45 bands, it would be incredibly complicated to display A, R, and T for each layer.  We need the averages Av-Rv-Tv and Ai-Ri-Ti to get a quick snapshot of what's going on, radiatively, with our different levels.  A good approach would be to display charts.  But I'm not getting into graphics here.  So here's a routine to find averages, no matter how many shortwave and longwave bands we have.  The code falls into the "tedious but straightforward" category:

    ' semigray finds two-tone radiative properties.
    Sub semigray()
        For L As Integer = 1 To ground
            With layer(L)
                .Av = 0.0#
                .Rv = 0.0#
                .Tv = 0.0#

                .Ai = 0.0#
                .Ri = 0.0#
                .Ti = 0.0#

                For b As Integer = swlo To swhi
                    .Av += .A(b)
                    .Rv += .R(b)
                    .Tv += .T(b)
                Next

                For b As Integer = lwlo To lwhi
                    .Ai += .A(b)
                    .Ri += .R(b)
                    .Ti += .T(b)
                Next

                .Av /= Nsw
                .Rv /= Nsw
                .Tv /= Nsw

                .Ai /= Nlw
                .Ri /= Nlw
                .Ti /= Nlw
            End With
        Next
    End Sub ' semigray

**36

	With semigray in place, we can now write a routine to show layer radiative properties.  And since we now have two display routines, let's call both from a general diagnostic routine:

    ' Main program.
    Sub Main()
        setup()
        diag()
    End Sub ' Main

    ' diag runs diagnostic routines.
    Sub diag()
        showPhys()
        showRad()
    End Sub ' diag

    ' showRad shows layer radiative properties.
    Sub showRad()
        Clear()
        WriteLine("L#     Av     Rv     Tv       Ai     Ri     Ti")
        WriteLine("--    -----  -----  -----    -----  -----  -----")
        semigray()

        For L As Integer = 1 To ground
            Write(L.ToString().PadLeft(2))

            With layer(L)
                Write(.Av.ToString("f3").PadLeft(9))
                Write(.Rv.ToString("f3").PadLeft(7))
                Write(.Tv.ToString("f3").PadLeft(7))

                Write(.Ai.ToString("f3").PadLeft(9))
                Write(.Ri.ToString("f3").PadLeft(7))
                Write(.Ti.ToString("f3").PadLeft(7))
            End With

            WriteLine()
        Next

        pause()
    End Sub ' showRad

**37

	Run what we've got now, and after the physical-values display, we see:

L#     Av     Rv     Tv       Ai     Ri     Ti
--    -----  -----  -----    -----  -----  -----
 1    0.000  0.003  0.997    0.150  0.000  0.850
 2    0.000  0.003  0.997    0.150  0.000  0.850
 3    0.000  0.003  0.997    0.150  0.000  0.850
 4    0.000  0.003  0.997    0.150  0.000  0.850
 5    0.000  0.003  0.997    0.150  0.000  0.850
 6    0.000  0.003  0.997    0.150  0.000  0.850
 7    0.000  0.003  0.997    0.150  0.000  0.850
 8    0.000  0.003  0.997    0.150  0.000  0.850
 9    0.000  0.003  0.997    0.150  0.000  0.850
10    0.000  0.003  0.997    0.150  0.000  0.850
11    0.000  0.003  0.997    0.150  0.000  0.850
12    0.000  0.003  0.997    0.150  0.000  0.850
13    0.000  0.003  0.997    0.150  0.000  0.850
14    0.000  0.003  0.997    0.150  0.000  0.850
15    0.000  0.003  0.997    0.150  0.000  0.850
16    0.000  0.003  0.997    0.150  0.000  0.850
17    0.000  0.003  0.997    0.150  0.000  0.850
18    0.000  0.003  0.997    0.150  0.000  0.850
19    0.000  0.003  0.997    0.150  0.000  0.850
20    0.000  0.003  0.997    0.150  0.000  0.850
21    0.860  0.140  0.000    0.996  0.004  0.000
Press [Enter] to go on.

	Don't forget to put calls to setRadAir and setRadGround in setup, or you'll see all zeros.

**38

	Now we need a way to display what happens to the layers as the model runs.  We also need to do the physics.  The first is easier, of course, so we'll tackle that first.  We can start off with some I/O routines.  VB's System.Console routines for locating the cursor and printing text are pretty verbose.  For a real number, it might be:

            SetCursorPosition(15, 3)
            Write(X.ToString("f3").PadLeft(7))

	I want something shorter, so I'll write three routines to put strings, integers, and real numbers where I want them on the screen.  I'll also change the inhuman start-with-0 Windows convention to the human-friendly start-with-1.
	All these routines have the same name, since VB allows overloading routines as long as they have differing "signatures"--order, number and type of passed arguments.
	Here r stands for row, c column, w field width, and d number of decimal places:

    ' put puts a string at row r, column c:
    Sub put(ByVal r As Integer, ByVal c As Integer, ByVal s As String)
        SetCursorPosition(c - 1, r - 1)
        Write(s)
    End Sub ' put

    ' Here's the integer version.
    Sub put(ByVal r As Integer, ByVal c As Integer, ByVal i As Integer, _
            ByVal w As Integer)
        SetCursorPosition(c - 1, r - 1)
        Write(i.ToString().PadLeft(w))
    End Sub ' put

    ' Here's the real-number version.
    Sub put(ByVal r As Integer, ByVal c As Integer, ByVal x As Double, _
            ByVal w As Integer, _
    ByVal d As Integer)
        Dim fmt As String = "f" & Trim(CStr(d))

        SetCursorPosition(c - 1, r - 1)
        Write(x.ToString(fmt).PadLeft(w))
    End Sub ' put

**39

	The main display screen will constantly change as we watch, until the simulation ends.  In addition to columns on the left showing what's happening in each layer, we will have 21 fields on the right describing overall results of the model.  Don't panic.  We won't put in most of those right away; they'll just display 0 until we gradually add them in as we go along.
	We will, however, need new global variables to hold the 0 values until we get to them.  Add:

    ' Earth
    Public Albedo As Double = 0.0#      ' Planetary albedo
    Public errPercent As Double = 0.0#  ' Energy conservation error, %
    Public F As Double = 0.0#           ' Absorbed climate flux density
    Public Fi As Double = 0.0#          ' Planetary IR radiated out
    Public Fv As Double = 0.0#          ' Planetary visual light reflected out
    Public H2Okg As Double = 0.0#       ' Water vapor burden
    Public sumF As Double = 0.0#        ' Total of Fi and Fv
    Public Te As Double = 0.0#          ' Planet emission temperature

    Public Tstrat As Double = 0.0#      ' Stratosphere temperature
    Public Ttrop As Double = 0.0#       ' Troposphere temperature
    Public TOAht As Double = 0.0#       ' Fictional Top-Of-Atmosphere height
    Public Tpause As Double = 0.0#      ' Tropopause height

    ' run
    Public cycle As Long = 0            ' How many times through main loop
    Public days As Double = 0.0#        ' Elapsed simulated time, days
    Public r2 As Double = 0.0#          ' Coefficient of determination
    Public res As Double = 99.0#        ' Temperature resolution achieved
    Public rms As Double = 0.0#         ' Root-mean-square error

**40

	Here's code to set the main screen up:

    ' setScreen prints labels for the main display.
    Sub setScreen()
        Const colA As Integer = 41

        Clear()
        WriteLine("L#  ht(m)  z (m)  lapse  T (K)")
        WriteLine("--  -----  -----  -----  -----")

        For L As Integer = 1 To ground
            put(L + 2, 1, L.ToString.PadLeft(2))
        Next

        put(1, colA, "Albedo              Fv")
        put(2, colA, "TOA                 Fi")
        put(2, 48, TOA, 5, 1)

        put(3, colA, "F                   sum")
        put(4, colA, "Te               K  error             %")

        put(6, colA, "cycle               res               K")
        put(7, colA, "days             d  time      '   " & Chr(34))
        put(8, colA, "H2O              kg")

        put(10, colA, "Tstrat           K  TOA ht            m")
        put(11, colA, "Ttrop            K  Tpause            m")
        put(12, colA, "Tsurf            K  ground        0   m")

        put(14, colA, "r2")
        put(15, colA, "rms              K")
    End Sub ' setScreen

**41

	Almost there.  Rewrite the main program:

    ' Main program.
    Sub Main()
        setup()
        diag()
        setScreen()

        cycle = 0
        days = 0.0#

        Do
            cycle += 1
            days += 1.0#

            transfer()
            heat()
            getRes()
            display()
        Loop Until res < 0.0001#

        pause()
    End Sub ' Main

	The algorithm should be clear.  We loop, calling four new routines repeatedly, until the temperature resolution falls to a very low level:  no layer has changed temperature by more than 1/10,000 of a degree.  transfer and heat do the physics, getRes is a simple check, and display shows the current values.  display is easiest, so we'll do that first.

**42

	The display routine is pretty straightforward:

    ' display shows the current state of the model.
    Sub display()
        Const colB As Integer = 48
        Const colD As Integer = 71

        For L = 1 To ground
            Dim row As Integer = L + 2
            put(row, 1, L, 2)

            With layer(L)
                put(row, 5, .height, 5, 0)
                put(row, 12, .z, 5, 0)
                put(row, 18, .gamma, 6, 2)
                put(row, 26, .temp, 5, 1)
            End With
        Next

        put(1, colB, Albedo, 5, 3)
        put(3, colB, F, 5, 1)
        put(4, colB, Te, 5, 1)

        put(1, colD, Fv, 5, 1)
        put(2, colD, Fi, 5, 1)
        put(3, colD, sumF, 5, 1)
        put(4, colD, errPercent, 5, 2)

        put(6, colB, cycle, 3)
        put(7, colB, days, 5, 1)

        SetCursorPosition(colB - 1, 8 - 1)
        Write(H2Okg.ToString("e2"))

        put(6, colD - 2, res, 7, 5)

        Dim mins As Integer = 0
        Dim secs As Integer = 0

        put(7, colD - 2, mins, 2)
        put(7, colD + 2, secs, 2)

        put(10, colB, Tstrat, 5, 1)
        put(11, colB, Ttrop, 5, 1)
        put(12, colB, layer(ground).temp, 5, 1)

        put(10, colD, TOAht, 5, 0)
        put(11, colD, Tpause, 5, 0)
    End Sub ' display

**43

	The transfer routine is the heart of the program.  It describes how radiative transfer works in the model, going from layer to layer.  It looks complicated, but in the end it's just bookkeeping according to the radiation identity--in each band, at each level, so much is absorbed, so much reflected, and the rest transmitted.  There's also emission, since the layers are warm.  Here's the code:

    ' transfer does radiative transfer
    Sub transfer()
        For b As Integer = 1 To Nbands
            Dim Fdown = band(b).solFrac * TOA
            Dim Fup = 0.0# ' Geothermal contribution is negligible.

            For L As Integer = 1 To ground
                With layer(L)
                    Dim fr As Double = D(band(b).hi, .temp) _
                    - D(band(b).lo, .temp)
                    .radOut(b) = fr * .A(b) * sb * .temp ^ 4

                    .above(b) = .A(b) * Fdown
                    .refUp(b) = .R(b) * Fdown
                    Fdown = .T(b) * Fdown + .radOut(b) + .refDown(b)
                End With
            Next L

            For L As Integer = ground To 1 Step -1
                With layer(L)
                    .below(b) = .A(b) * Fup
                    .refDown(b) = .R(b) * Fup
                    Fup = .T(b) * Fup + .radOut(b) + .refUp(b)
                End With
            Next L
        Next b
    End Sub ' transfer

	b loops through bands, L loops through layers.

**44

	For heat, we need an equation for heat content.  Heat, the most general form of energy, is measured in joules in the SI.  An object's heat content is

	H = m cp T										|29

where m is the object's mass (in kilograms in the SI), cp the specific heat capacity (J kg-1 K-1), T the temperature (K).  Now you know why I added the seemingly redundant variable cpm in the layer structure.  No point repeating the same calculation with the same values every time through the loop.
	Since the area of each layer doesn't increase much with altitude--the altitude isn't great compared with Earth's radius--we can take each layer as being roughly the same surface area as the Earth.  Given that, and the equation above, here's the code for the heating routine.  Again, it's pretty much just energy bookkeeping--so much absorbed by each layer in each band, so much emitted, net change causes temperature change:

    ' heat sets the temperature of each layer.
    Sub heat()
        Const factor As Double = 86400.0# * area

        For L As Integer = 1 To ground
            Dim HeatIn As Double = 0.0#
            Dim HeatOut As Double = 0.0#

            With layer(L)
                For b As Integer = 1 To Nbands
                    HeatIn += .above(b) + .below(b)

                    If L = ground Then
                        HeatOut += .radOut(b) ' Up only.
                    Else
                        HeatOut += 2.0# * .radOut(b) ' Up and down.
                    End If
                Next

                Dim net As Double = (HeatIn - HeatOut) * factor
                .oldTemp = .temp
                .temp += (net / .cpm)
            End With
        Next
    End Sub ' heat

	Where did "86400.0#" come from?  Flux density is in watts per square meter, which is joules per square meter per second.  We're using a one-day time step.  There are 60 x 60 x 24 = 86,400 seconds in a day.  And area square meters in a layer.  So the total energy the layer absorbs in one day is the net flux density, times the area, times the elapsed time in seconds.

**45

	The last unfinished routine, getRes, is trivial.  It just finds the biggest temperature change in any layer:

    ' getRes finds the temperature resolution.
    Sub getRes()
        res = 0.0# ' Assume everything's perfect.

        For L As Integer = 1 To ground
            Dim delta As Double = Abs(layer(L).temp - layer(L).oldTemp)

            If delta > res Then
                res = delta
            End If
        Next
    End Sub ' getRes

**46

	Run the program and we get a really cool display, with the temperatures in each layer changing as we watch.  It settles down after a minute, when the resolution reaches the critical level, and we stop.

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.000        Fv          0.0
--  -----  -----  -----  -----          TOA    340.3        Fi          0.0
 1      0      0   0.00  227.1          F        0.0        sum         0.0
 2      0      0   0.00  234.9          Te       0.0     K  error      0.00   %
 3      0      0   0.00  242.1
 4      0      0   0.00  248.7          cycle  357          res     0.00010   K
 5      0      0   0.00  254.8          days   357.0     d  time     0'  0"
 6      0      0   0.00  260.5          H2O    0.00e+000 kg
 7      0      0   0.00  265.8
 8      0      0   0.00  270.9          Tstrat   0.0     K  TOA ht        0   m
 9      0      0   0.00  275.7          Ttrop    0.0     K  Tpause        0   m
10      0      0   0.00  280.2          Tsurf  335.9     K  ground        0   m
11      0      0   0.00  284.5
12      0      0   0.00  288.7          r2
13      0      0   0.00  292.7          rms              K
14      0      0   0.00  296.5
15      0      0   0.00  300.2
16      0      0   0.00  303.7
17      0      0   0.00  307.2
18      0      0   0.00  310.5
19      0      0   0.00  313.8
20      0      0   0.00  316.9
21      1     -1   0.00  335.9

	After this I'll leave off the final "Press [Enter] to go on."
	The Earth is warm--336 K at the surface, 48 K higher than it should be.  This is partly from low albedo (no clouds), and partly from the lack of convection (we're tracing radiation only).  But before we improve the model, let’s improve the diagnostics.

**47

	Add a new global variable to the "Earth" section:

    Public toSpace(Nbands) As Double    ' Emitted in each band.

	And use it in transfer, just after the second layer loop and inside the band loop:

            toSpace(b) = Fup

	Then we can write a routine to check conservation of energy:

    ' conserve checks energy conservation.
    Sub conserve()
        Fv = 0.0#
        Fi = 0.0#

        For b As Integer = swlo To swhi
            Fv += toSpace(b)
        Next

        For b As Integer = lwlo To lwhi
            Fi += toSpace(b)
        Next

        sumF = Fi + Fv
        Albedo = Fv / TOA
        F = TOA * (1.0# - Albedo)
        Te = (F / sb) ^ 0.25

        errPercent = 100.0# * sumF / TOA - 100.0#
    End Sub ' conserve

**48

	Call conserve from within the display routine, right after the big layer loop.  Run the program again.  When it ends, the right-hand section reads:

                                        Albedo 0.193        Fv         65.5
                                        TOA    340.3        Fi        274.6
                                        F      274.7        sum       340.1
                                        Te     263.8     K  error     -0.044  %

                                        cycle  357          res     0.00010   K
                                        days   357.0     d  time     0'  0"
                                        H2O    0.00e+000 kg

                                        Tmeso    0.0     K  TOA ht        0   m
                                        Ttrop    0.0     K  Tpause        0   m
                                        Tsurf  335.9     K  ground        0   m

                                        r2
                                        rms              K

	Our albedo is indeed low, due to the lack of clouds.  But we're conserving energy to within 0.044%, which is a good sign.
	This is a very crude model indeed.  But don't knock it, we've achieved something--a real, honest-to-God time-marching radiative equilibrium model of Earth's climate system.  It's not accurate, it leaves out all kinds of things, but if Earth had no clouds and if its atmosphere were semigray, and if there were no convection or conduction, this is what it would be like.
	Now to make things gradually more realistic.

**49

	First, let's make the display more readable.  It will help to fill in the column values.  Here's a routine to find heights and lapse rates:

    ' lapse finds lapse rates and heights.
    Sub lapse()
        For L As Integer = 1 To N
            With layer(L)
                Dim V As Double = .num * Rstar * .temp / .p
                .height = V / area
            End With
        Next

        Dim baseAltitude As Double = 0.0#

        For L As Integer = N To 1 Step -1
            With layer(L)
                .z = 0.5# * .height + baseAltitude
                baseAltitude += .height

                Dim M As Integer = L + 1
                Dim dT As Double = .temp - layer(M).temp
                Dim dz As Double = .z - layer(M).z
                .gamma = -1000.0# * dT / dz
            End With
        Next
    End Sub ' lapse

**50

	Call it from the main routine just before display.  Then we run the program and get:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.193        Fv         65.5
--  -----  -----  -----  -----          TOA    340.3        Fi        274.6
 1  13304  29354   0.88  227.1          F      274.7        sum       340.1
 2   4589  20408   1.93  234.9          Te     263.8     K  error     -0.044  %
 3   2837  16695   2.68  242.1
 4   2082  14235   3.26  248.7          cycle  357          res     0.00010   K
 5   1659  12365   3.74  254.8          days   357.0     d  time     0'  0"
 6   1388  10842   4.13  260.5          H2O    0.00e+000 kg
 7   1198   9549   4.47  265.8
 8   1058   8421   4.76  270.9          Tmeso    0.0     K  TOA ht        0   m
 9    950   7417   5.01  275.7          Ttrop    0.0     K  Tpause        0   m
10    864   6510   5.23  280.2          Tsurf  335.9     K  ground        0   m
11    794   5681   5.43  284.5
12    735   4916   5.60  288.7          r2
13    686   4205   5.76  292.7          rms              K
14    643   3541   5.90  296.5
15    607   2916   6.03  300.2
16    574   2325   6.15  303.7
17    545   1766   6.26  307.2
18    520   1233   6.36  310.5
19    497    725   6.45  313.8
20    476    238  79.88  316.9
21      1     -1   0.00  335.9

	The lapse rate is defined as the negative of the temperature change with altitude:

	Γ = -dT / dz									|30

	Since the convention is to measure it in K km-1, I multiply by 1,000 in the lapse rate routine.  Otherwise we'd get K m-1.
	Note the very high lapse rate, 80 K/km, in the lowest layer of the troposphere.  This is unrealistically high.  In reality, a lapse rate this high would cause convective overturning, and heat would mix between layers.  But we're not doing convection yet, so we have to tolerate this for now.

**51

	Here's a routine to fill in data for the meteorological layers.

    ' metLayers finds meteorological layer data.
    Sub metLayers()
        Const topStrat As Integer = 1
        Const bottomStrat As Integer = 4
        Const Nstrat As Integer = bottomStrat - topStrat + 1

        Const topTrop As Integer = bottomStrat + 1
        Const bottomTrop As Integer = N
        Const Ntrop As Integer = bottomTrop - topTrop + 1

        Tstrat = 0.0#
        Ttrop = 0.0#

        For L As Integer = topStrat To bottomStrat
            Tstrat += layer(L).temp
        Next

        For L As Integer = topTrop To bottomTrop
            Ttrop += layer(L).temp
        Next

        Tstrat /= Nstrat
        Ttrop /= Ntrop

        TOAht = layer(topStrat).z + 0.5# * layer(topStrat).height
        Tpause = layer(topTrop).z + 0.5# * layer(topTrop).height
    End Sub ' metLayers

**52

	In display, call metLayers right after the call to conserve.  Now the right-hand section reads:

                                        Albedo 0.193        Fv         65.5
                                        TOA    340.3        Fi        274.6
                                        F      274.7        sum       340.1
                                        Te     263.8     K  error     -0.044  %

                                        cycle  357          res     0.00010   K
                                        days   357.0     d  time     0'  0"
                                        H2O    0.00e+000 kg

                                        Tstrat 238.2     K  TOA ht    36006   m
                                        Ttrop  288.9     K  Tpause    13194   m
                                        Tsurf  335.9     K  ground        0   m

                                        r2
                                        rms              K

	All the layers are too hot.  The tropopause, due to the hot troposphere, is higher than it should be (11 km on average), though not by much.  But the top-of-atmosphere is way too low.  This is an unavoidable consequence of quantizing the atmosphere into distinct layers.  If we increased the number of layers, we could do better, but for now, let's ignore this.  We're using pressure as our vertical coordinate.

**53

	The next step is to set each air layer's absorptance, reflectance, and transmittance from real physical principles, based on the gases it contains.  This will be a long, complicated process, I'm afraid.  But we can split it into two major parts.  First, distributing the gases among the layers.  Second, setting the radiative properties.  Which means we'll need a real band scheme, not just the shortwave and longwave averages we've used so far.

**54

	To parcel out gases, here is the code to the distrib routine.  It isn't complicated:

    ' distrib parcels out gases to different layers.
    Sub distrib()
        For L As Integer = 1 To N
            With layer(L)
                For gs As Integer = 1 To Ngases
                    Select Case gs
                        Case H2O : .SM(H2O) = water(L)
                        Case O3 : .SM(O3) = ozone(L)
                        Case Else : .SM(gs) = absorber(gs).q * .mass
                    End Select

                    .SM(gs) /= area
                Next
            End With
        Next
    End Sub ' distrib

	Of course, this requires two more functions, to distribute ozone and water vapor, respectively.  For the moment, we can make these "stub" functions without meaningful content:

    ' ozone distributes O3 to different layers.
    Function ozone(ByVal L As Integer) As Double
        Return 0.0#
    End Function ' ozone

    ' water distributes H2O to different layers.
    Function water(ByVal L As Integer) As Double
        Return 0.0#
    End Function ' water

**55

	Ozone has a peculiar distribution.  Most of it is in the stratosphere--it gets formed there by sunlight breaking apart oxygen molecules into oxygen atoms:

	O2 + ν => O + O									|31
	O2 + O => O3										|32

	There are also reactions, involving third-party elements, that break up ozone into oxygen again.  The net reaction is:

	O3 + ν => O2 + O									|33

	The balance between the two types of reaction determines the equilibrium level of ozone.  About 90% of atmospheric ozone is in the stratosphere, and the rest is in the troposphere, a little thicker near the ground due to pollution-related photochemistry.

**56

	The U.S. Standard Atmosphere put together in 1976 by NOAA, NASA, and the US Air Force is a standard for comparing atmosphere models.  One useful thing it does is provide ozone fractions at each altitude.  Each altitude has an associated pressure level.  By matching the pressure levels, I put together a table of ozone mass in each layer.  Thus the "ozone" function:

    ' ozone returns the O3 mass in a layer.
    Function ozone(ByVal L As Integer) As Double
        Dim O3mass() As Double = {0.0#, _
        1879000000000.0#, 701600000000.0#, _
        266300000000.0#, 153800000000.0#, _
        93040000000.0#, 51370000000.0#, _
        33440000000.0#, 24200000000.0#, _
        20570000000.0#, 17800000000.0#, _
        16680000000.0#, 15540000000.0#, _
        14780000000.0#, 14600000000.0#, _
        14410000000.0#, 14220000000.0#, _
        13620000000.0#, 12990000000.0#, _
        12420000000.0#, 11900000000.0#}

        Return O3mass(L)
    End Function ' ozone

	The total mass of ozone in the atmosphere, from this analysis, is about 3.38 trillion kilograms.  Walker (1977) estimated 3.3 trillion, and modern values may be closer to 3.0 trillion due to ozone depletion.  Whatever.  We're going to try to match our models against the USSA-76, so we'll stick with this distribution for now.

**57

	Since temperatures at Earth's surface and in the troposphere cycle around the freezing point of water, water is a volatile in Earth's atmosphere--easily changing phase, solid to liquid to gas or whatever.  Water, mostly from the oceans, evaporates, condenses as it rises, precipitates out as rain, snow, etc.  Und so weiter.  World without end, amen.
	The water vapor present is strongly related to the ambient temperature by the Clausius-Clapeyron law, discovered and explained in the 19th century:

	dP    L
	 = 										|34
	dT   T dV

where P is the "saturation vapor pressure," T temperature, L the latent heat for the transition from liquid to gas, and V the "specific volume," or reciprocal of the density.  d means "small change in," and of course this is really a differential equation.  And in practice, it's an extremely difficult equation to work with.  Meteorologists, climatologists and chemists usually rely on approximations.  Of course, approximations can be made as accurate as needed.

**58

	A good approximation for the Clausius-Clapeyron relation is the August-Roche-Magnus formula:

	                  17.625 Tc
	eS = 6.1094 exp []						|35
	                 Tc + 243.04

where saturation vapor pressure eS is in millibars or "hektopascals" (hPa), and temperature Tc is in degrees Celsius.  Since we're dealing with pressure in pascals and temperature in Kelvin, we make some algebraic substitutions:

	                17.625 T - 4814.26875
	es = 610.94 exp[]					|36
	                       T - 30.11

which gives us a function for saturation vapor pressure:

    ' cc approximates the Clausius-Clapeyron relation
    ' with the August-Roche Magnus formula.
    Function cc(ByVal L As Integer) As Double
        Const k0 As Double = 610.94#
        Const k1 As Double = 17.625#
        Const k2 As Double = 4814.26875#
        Const k3 As Double = 30.11#

        Dim Tk As Double = layer(L).temp
        Dim es As Double = k0 * Exp((k1 * Tk - k2) / (Tk - k3))

        Return es
    End Function ' cc

**59

	Now we can put together a simple function for the water vapor mass in a given level:

    ' water distributes H2O to different layers.
    Function water(ByVal L As Integer) As Double
        With layer(L)
            Dim pH2O As Double = .RH * cc(L)
            Dim XH2O As Double = pH2O / .p
            Dim qH2O As Double = 0.622# * XH2O
            Return qH2O * .mass
        End With
    End Function ' water

	The actual pressure is the relative humidity times the saturation pressure.  The constant 0.622 is the ratio of the molecular weight of water vapor to that of air (taking the mean air value as 28.94), so multiplying volume fraction X by that gives mass fraction q.

**60

	Let's call distrib from setup, right after the call to setPhysGround.  In addition, let's add another diagnostic display, so we can see what we're getting:

    ' showSM shows mass paths.
    Sub showSM()
        Clear()
        WriteLine("L#    CO2      H2O       O3")
        WriteLine("--  -------  -------  ---------")

        For L As Integer = 1 To N
            Write(L.ToString.PadLeft(2))

            With layer(L)
                Write(.SM(CO2).ToString("f3").PadLeft(9))
                Write(.SM(H2O).ToString("f3").PadLeft(9))
                Write(.SM(O3).ToString("f5").PadLeft(11))
            End With

            WriteLine()
        Next

        pause()
    End Sub ' showSM

	Call that from diag right after the call to ShowPhys.  The second screen now shows us:

L#    CO2      H2O       O3
--  -------  -------  ---------
 1    0.247    0.067    0.00368
 2    0.247    0.271    0.00138
 3    0.247    0.312    0.00052
 4    0.247    0.329    0.00030
 5    0.247    0.339    0.00018
 6    0.247    0.345    0.00010
 7    0.247    0.349    0.00007
 8    0.247    0.352    0.00005
 9    0.247    0.355    0.00004
10    0.247    0.356    0.00003
11    0.247    0.358    0.00003
12    0.247    0.359    0.00003
13    0.247    0.360    0.00003
14    0.247    0.361    0.00003
15    0.247    0.362    0.00003
16    0.247    0.363    0.00003
17    0.247    0.363    0.00003
18    0.247    0.364    0.00003
19    0.247    0.364    0.00002
20    0.247    0.365    0.00002

	If you get all zeros in the CO2 column, you forgot to call setAbsorbers from the master setup routine.

**61

	Okay.  We've got a more realistic atmosphere.  Now we need a serious band structure.
	Before we go any further, let's add yet another diagnostic screen.  This one will show the fraction of solar energy in each band, plus the fraction of energy from a blackbody at 255 K, typical of the atmosphere:

    ' showSolar shows solar energy distribution.
    Sub showSolar()
        Clear()
        WriteLine("b#   low    high  solar    f255")
        WriteLine("--  -----  -----  ------  ------")

        For b As Integer = 1 To Nbands
            Write(b.ToString.PadLeft(2))

            With band(b)
                Write(.lo.ToString("f1").PadLeft(7))
                Write(.hi.ToString("f1").PadLeft(7))
                Write(.solFrac.ToString("f4").PadLeft(8))

                Dim f255 As Double = D(.hi, 255.0#) - D(.lo, 255.0#)

                Write(f255.ToString("f4").PadLeft(8))
            End With

            WriteLine()
        Next

        WriteLine()
        Write("Press [Enter] to go on.  ")
        Dim st As String = ReadLine()
    End Sub ' showSolar

	This new third screen now shows:

b#   low    high  solar    f255
--  -----  -----  ------  ------
 1    0.2    0.2  0.0039  0.0000
 2    0.2    0.2  0.0041  0.0000
 3    0.2    0.3  0.0126  0.0000
 4    0.3    0.3  0.0078  0.0000
 5    0.3    0.3  0.0093  0.0000
 6    0.3    0.3  0.0070  0.0000
 7    0.3    0.4  0.0767  0.0000
 8    0.4    0.7  0.3665  0.0000
 9    0.7    1.8  0.4334  0.0000
10    1.8    1.9  0.0096  0.0000
11    1.9    2.0  0.0081  0.0000
12    2.0    2.1  0.0069  0.0000
13    2.1    2.2  0.0059  0.0000
14    2.2    2.5  0.0133  0.0000
15    2.5    2.6  0.0033  0.0000
16    2.6    2.9  0.0077  0.0000
17    2.9    3.0  0.0020  0.0000
18    3.0    3.5  0.0072  0.0001
19    3.5    4.1  0.0049  0.0005
20    4.1    4.5  0.0021  0.0009
21    4.5    5.0  0.0018  0.0023
22    5.0    6.5  0.0027  0.0209
23    6.5    7.0  0.0005  0.0130
24    7.0    8.0  0.0006  0.0355
25    8.0    9.1  0.0004  0.0498
26    9.1    9.3  0.0001  0.0097
27    9.3    9.3  0.0000  0.0029
28    9.3    9.4  0.0000  0.0042
29    9.4    9.5  0.0000  0.0040
30    9.5    9.6  0.0000  0.0072
31    9.6    9.7  0.0000  0.0073
32    9.7    9.9  0.0000  0.0060
33    9.9   10.0  0.0000  0.0052
34   10.0   10.0  0.0000  0.0053
35   10.0   10.1  0.0000  0.0050
36   10.1   12.0  0.0003  0.1071
37   12.0   13.2  0.0001  0.0696
38   13.2   13.5  0.0000  0.0169
39   13.5   14.0  0.0000  0.0242
40   14.0   16.0  0.0001  0.0978
41   16.0   16.7  0.0000  0.0289
42   16.7   25.0  0.0001  0.2432
43   25.0   30.0  0.0000  0.0740
44   30.0   40.0  0.0000  0.0768
45   40.0  125.0  0.0000  0.0779

	Note that 99% of solar energy is at 4 μ or less.  Note also that only 10% of atmospheric energy, roughly, falls in the CO2 15-μ band.  Next, we implement the band scheme.

**62

	Our first step is to add five (5, count 'em, 5!) new arrays to the absorber structure:

        Dim g() As Double   ' Asymmetry factor.
        Dim ka() As Double  ' Absorption coefficient.
        Dim ke() As Double  ' Extinction coefficient.
        Dim ks() As Double  ' Scattering coefficient.
        Dim om() As Double  ' Single-scattering albedo.

	Along with two real numbers:

        Dim m As Double     ' Pressure exponent.
        Dim n As Double     ' Temperature exponent.

	g, the asymmetry factor, relates how much light is scattered forward as compared to backward, and ranges from -1 to 1.  ka, ke, and ks are just variants on the absorption coefficient we noted when we discussed Beer's Law.  Extinction is the overall effect, and it can happen in two ways:  light can be absorbed, or it can be scattered--in a one-dimensional setup, reflected.  om is actually omega (ω), the "single-scattering albedo," which is the fraction of extinction accounted for by scattering:

	ω = ks / ke										|37

	But again, don't panic--we're not going to use most of these until much later, when we get to clouds.  The only one we care about for the moment is ka, the absorption coefficient.
	m and n are exponents, used to factor in the effect of layer pressure and temperature on band strength.  We'll explore their use later.

**63

	Add more code to dimension:

        For a As Integer = 1 To Nabs
            With absorber(a)
                ReDim .g(Nbands)
                ReDim .ka(Nbands)
                ReDim .ke(Nbands)
                ReDim .ks(Nbands)
                ReDim .om(Nbands)
            End With
        Next

	And at the end of setBands, add yet more code:

        For b As Integer = 1 To Nbands
            absorber(CH4).ka(b) = 0.0#
            absorber(CO2).ka(b) = 0.0#
            absorber(H2O).ka(b) = 0.0#
            absorber(O3).ka(b) = 0.0#
        Next b

	For the moment, we will let all the absorption coefficients be zero, so we see what happens with no greenhouse effect.

**64

	This is the biggest change we need to make.  We completely rewrite the big loop in setRadAir that finds absorptance, reflectance, and transmittance for each layer in each band:

    ' setRadAir sets air layer radiative properties.
    Sub setRadAir()
        For L = 1 To N
            With layer(L)
                For b As Integer = 1 To Nbands
                    Dim Tch4 As Double = Exp(-tau(CH4, L, b))
                    Dim Tco2 As Double = Exp(-tau(CO2, L, b))
                    Dim Th2o As Double = Exp(-tau(H2O, L, b))
                    Dim To3 As Double = Exp(-tau(H2O, L, b))

                    Dim Tclear As Double = Tch4 * Tco2 * Th2o * To3

                    If b <= swhi Then
                        .R(b) = 0.06# * .mfrac
                    Else
                        .R(b) = 0.0#
                    End If

                    If (Tclear + .R(b)) >= 1.0 Then
                        .A(b) = 0.0#
                        .T(b) = 1.0# - .R(b)
                    Else
                        .T(b) = Tclear
                        .A(b) = 1.0# - .R(b) - .T(b)
                    End If
                Next b
            End With
        Next L
    End Sub ' setRadAir

	We will need a tau routine to find optical thickness due to each absorber.  We'll do that later.  For the moment, note that I am using the "multiplicative approximation"--that net transmittance is the product of all the individual transmittances.  This is not always accurate, and the literature is filled with different mathematical approaches to "band overlapping."  But it's a good approximation in most cases.

**65

	Here's the optical depth routine:

    ' tau finds optical depth due to an absorber.
    Function tau(ByVal a As Integer, ByVal L As Integer, ByVal b As Integer) As Double
        Const diff = 1.66#

        Dim pFactor As Double = 0.0#
        Dim tFactor As Double = 0.0#
        Dim uStar As Double = 0.0#

        With layer(L)
            pFactor = (.p / P0) ^ absorber(a).m
            tFactor = (T0 / .temp) ^ absorber(a).n
            uStar = .SM(a) * pFactor * tFactor
        End With

        Return diff * absorber(a).ka(b) * uStar
    End Function ' tau

	Allow me to explain.
	Absorption lines can be broadened by several mechanisms.  In the lower atmosphere (most of the mass), "pressure broadening" is the major mechanism.  Higher up, it's "Doppler broadening" due to molecular motion.  The result is that greenhouse gases are more effective at high pressure and low temperature.  Atmosphere physicists express it this way:

	u* = (P /P0)m (T0 / T)n								|38

The symbol u is often used for specific mass or mass path (our SM in the layerType structure).  u* is then the "scaled mass path," as defined by the above equation.  P and T are pressure and temperature.  P0 and T0 are reference values of pressure and temperature, usually taken as the sea-level mean values of 101,325 Pa and 288.15 K, respectively.  m and n are exponents.  Many researchers find different exponents for each greenhouse gas, but I'm using a primitive model here (Hart 1978) where m is 0.5 and n is 0.25.  A lot of models use n = 1/2 m.
	To go on.  1.66 is the "diffusion coefficient."  Optical depth χ for a three-dimensional setup can be translated to τ for a 1-D setup by simply multiplying by 5/3.  I know, 1.67 would be a better approximation, but for some reason a lot of climatologists use 1.66.  Go figure.

**66

	We run.  The two-tone radiative properties are:

L#     Av     Rv     Tv       Ai     Ri     Ti
--    -----  -----  -----    -----  -----  -----
 1    0.000  0.003  0.997    0.000  0.000  1.000
 2    0.000  0.003  0.997    0.000  0.000  1.000
 3    0.000  0.003  0.997    0.000  0.000  1.000
 4    0.000  0.003  0.997    0.000  0.000  1.000
 5    0.000  0.003  0.997    0.000  0.000  1.000
 6    0.000  0.003  0.997    0.000  0.000  1.000
 7    0.000  0.003  0.997    0.000  0.000  1.000
 8    0.000  0.003  0.997    0.000  0.000  1.000
 9    0.000  0.003  0.997    0.000  0.000  1.000
10    0.000  0.003  0.997    0.000  0.000  1.000
11    0.000  0.003  0.997    0.000  0.000  1.000
12    0.000  0.003  0.997    0.000  0.000  1.000
13    0.000  0.003  0.997    0.000  0.000  1.000
14    0.000  0.003  0.997    0.000  0.000  1.000
15    0.000  0.003  0.997    0.000  0.000  1.000
16    0.000  0.003  0.997    0.000  0.000  1.000
17    0.000  0.003  0.997    0.000  0.000  1.000
18    0.000  0.003  0.997    0.000  0.000  1.000
19    0.000  0.003  0.997    0.000  0.000  1.000
20    0.000  0.003  0.997    0.000  0.000  1.000
21    0.860  0.140  0.000    0.996  0.004  0.000

	The final result is:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.181        Fv         61.6
--  -----  -----  -----  -----          TOA    340.3        Fi        278.5
 1  14941  29579   0.00  255.0          F      278.6        sum       340.1
 2   4980  19618   0.00  255.0          Te     264.8     K  error     -0.041  %
 3   2988  15634   0.00  255.0
 4   2134  13072   0.00  255.0          cycle   98          res     0.00009   K
 5   1660  11175   0.00  255.0          days    98.0     d  time     0'  0"
 6   1358   9666   0.00  255.0          H2O    0.00e+000 kg
 7   1149   8412   0.00  255.0
 8    996   7339   0.00  255.0          Tstrat 255.0     K  TOA ht    37049   m
 9    879   6402   0.00  255.0          Ttrop  255.0     K  Tpause    12005   m
10    786   5569   0.00  255.0          Tsurf  265.3     K  ground        0   m
11    711   4820   0.00  255.0
12    650   4140   0.00  255.0          r2
13    598   3516   0.00  255.0          rms              K
14    553   2940   0.00  255.0
15    515   2406   0.00  255.0
16    482   1908   0.00  255.0
17    453   1440   0.00  255.0
18    427   1000   0.00  255.0
19    404    585   0.00  255.0
20    383    192  53.53  255.0
21      1     -1   0.00  265.3

	As expected:  Very cold.  In fact, the Earth, without greenhouse gases, is frozen over.  As Fourier pointed out in 1824.

**67

	Let's add just one absorption band to the atmosphere:  the 15-micron band of CO2.  At the end of setBands, we call a new routine:

        setkCO2()

which is defined thusly:

    ' setkCO2 sets ka for carbon dioxide.
    Sub setkCO2()
        With absorber(CO2)
            .ka(40) = 20.0#
        End With
    End Sub ' setkCO2

	Our radiative screen shows:

L#     Av     Rv     Tv       Ai     Ri     Ti
--    -----  -----  -----    -----  -----  -----
 1    0.000  0.003  0.997    0.028  0.000  0.972
 2    0.000  0.003  0.997    0.035  0.000  0.965
 3    0.000  0.003  0.997    0.036  0.000  0.964
 4    0.000  0.003  0.997    0.037  0.000  0.963
 5    0.000  0.003  0.997    0.038  0.000  0.962
 6    0.000  0.003  0.997    0.038  0.000  0.962
 7    0.000  0.003  0.997    0.038  0.000  0.962
 8    0.000  0.003  0.997    0.038  0.000  0.962
 9    0.000  0.003  0.997    0.038  0.000  0.962
10    0.000  0.003  0.997    0.038  0.000  0.962
11    0.000  0.003  0.997    0.038  0.000  0.962
12    0.000  0.003  0.997    0.038  0.000  0.962
13    0.000  0.003  0.997    0.038  0.000  0.962
14    0.000  0.003  0.997    0.038  0.000  0.962
15    0.000  0.003  0.997    0.038  0.000  0.962
16    0.000  0.003  0.997    0.038  0.000  0.962
17    0.000  0.003  0.997    0.038  0.000  0.962
18    0.000  0.003  0.997    0.038  0.000  0.962
19    0.000  0.003  0.997    0.038  0.000  0.962
20    0.000  0.003  0.997    0.038  0.000  0.962
21    0.860  0.140  0.000    0.996  0.004  0.000

	It takes a long time to run--sorry about that.  Finally, when it settles down, we see:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.181        Fv         61.7
--  -----  -----  -----  -----          TOA    340.3        Fi        278.4
 1   8379  21253   2.55  143.0          F      278.5        sum       340.1
 2   3078  15525   4.95  157.6          Te     264.7     K  error     -0.032  %
 3   1994  12989   5.98  170.1
 4   1512  11237   6.56  180.6          cycle  5874         res     0.00010   K
 5   1234   9863   6.97  189.6          days   5874.0    d  time     0'  0"
 6   1052   8720   7.27  197.6          H2O    0.00e+000 kg
 7    923   7732   7.53  204.8
 8    826   6858   7.75  211.4          Tstrat 162.8     K  TOA ht    25442   m
 9    749   6070   7.94  217.5          Ttrop  233.5     K  Tpause    10481   m
10    688   5352   8.12  223.2          Tsurf  271.7     K  ground        0   m
11    638   4689   8.28  228.5
12    595   4072   8.43  233.6          r2
13    559   3495   8.58  238.5          rms              K
14    528   2952   8.72  243.2
15    500   2438   8.85  247.7
16    476   1949   8.98  252.0
17    455   1484   9.10  256.2
18    436   1039   9.22  260.2
19    418    612   9.33  264.1
20    403    201  18.58  268.0
21      1     -1   0.00  271.7

	CO2 has raised the surface temperature from 265.3 K to 271.7 K, a 6.4 K increase.  Not enough to unfreeze the ground, but substantial.  	Sure did take a long time to run, though.  Five minutes?  Ten?  Let's start using those time fields.

**68

	Add a new global variable:

    Public startTime As DateTime        ' When simulation starts.

	Just before the big loop in Main, add:

        startTime = Now

	Finally, in display, just after the call to metLayers, add:

        Dim seconds As Double = Now.Subtract(startTime).TotalSeconds
        Dim mins As Integer = CInt(seconds) \ 60
        Dim secs As Integer = CInt(seconds) Mod 60

	Delete the declarations put down earlier, below this, for mins and secs.
	Now, when we run the program, we see elapsed time counting off second by second.  This version, with one absorption band, takes 13 minutes and 13 seconds to run.  Yikes!
	Next, let's go to a full-scale, professional-style band scheme.  No more putzing around with no greenhouse gases or single absorption bands.

**69

	Here's the new band scheme.  The first table shows the shortwave scheme, the second one shows longwave.

low
	high
	kCH4
	kCO2
	kH2O
	kO3

0.175	0.225				1527.0	
0.225	0.245				9384.1	
0.245	0.280				15132.0	
0.280	0.295				2147.0	
0.295	0.310				355.3	
0.310 	0.320				62.6	
0.320	0.400				1.50	
0.4	0.7				2.51	
0.7	1.8					
1.8	1.9			714		
1.9	2.0		334	714		
2.0	2.1		334			
2.1	2.2					
2.2	2.5	147				
2.5	2.6			71.4		
2.6	2.9		71	71.4		
2.9	3.0			71.4		
3.0	3.5					
3.5
	4.1
	18

low
	high
	kCH4
	kCO2
	kH2O
	kO3

4.1	4.5	18.0	9.356			
4.5	5.0					
5.0	6.5			11		
6.5	7.0	3.2		11		
7.0	8.0			11		
8.000	9.075					
9.075	9.265				34.0	
9.265	9.320				130.0	
9.320	9.400				272.0	
9.400	9.475				406.0	
9.475	9.610				395.0	
9.610	9.745				333.0	
9.745	9.855				216.0	
9.855	9.950				112.0	
9.950	10.045				40.0	
10.045	10.135				7.5	
10.135	12.000					
12.000	13.235			0.86		
13.235	13.545		0.17	0.86		
13.545	14.00		0.82	0.86		
14.00	16.00		20.00	0.86		
16.00	16.66		0.88	0.86		
16.66	25.00			0.86		
25.00	30.00			4.00		
30.00	40.00			16.00		
40.00
	125.00
			40.00

	Methane is from Essenhigh (2001), carbon dioxide from Essenhigh and Gonima (1992), water vapor from Essenhigh and Hanel et al. (1963), and ozone from Gonima and Chou and Lee (1996).

**70

	Here's new code to put at the end of setBands:

        setkCH4()
        setkCO2()
        setkH2O()
        setkO3()

	And here are those subroutines:

    ' setkCH4 sets ka for methane.
    Sub setkCH4()
        With absorber(CH4)
            .ka(14) = 147.0#
            .ka(19) = 18.0#
            .ka(20) = 18.0#
            .ka(23) = 3.2#
        End With
    End Sub ' setkCH4

    ' setkCO2 sets ka for carbon dioxide.
    Sub setkCO2()
        With absorber(CO2)
            .ka(11) = 334.0#
            .ka(12) = 334.0#
            .ka(16) = 71.0#
            .ka(20) = 9.356#

            .ka(38) = 0.17#
            .ka(39) = 0.82#
            .ka(40) = 20.0#
            .ka(41) = 0.88#
        End With
    End Sub ' setkCO2

    ' setkH2O sets ka for water vapor.
    Sub setkH2O()
        With absorber(H2O)
            .ka(10) = 714.0#
            .ka(11) = 714.0#
            .ka(15) = 71.4#
            .ka(16) = 71.4#
            .ka(17) = 71.4#
            .ka(22) = 11.0#
            .ka(23) = 11.0#
            .ka(24) = 11.0#

            For b As Integer = 37 To 42
                .ka(b) = 0.86#
            Next

            .ka(43) = 4.0#
            .ka(44) = 16.0#
            .ka(45) = 40.0#
        End With
    End Sub ' setkH2O

    ' setkO3 sets ka for ozone.
    Sub setkO3()
        With absorber(O3)
            .ka(1) = 1527.0#
            .ka(2) = 9384.1#
            .ka(3) = 15132.0#
            .ka(4) = 2147.0#
            .ka(5) = 355.3#
            .ka(6) = 62.6#
            .ka(7) = 1.5#
            .ka(8) = 2.51#

            .ka(26) = 34.0#
            .ka(27) = 130.0#
            .ka(28) = 272.0#
            .ka(29) = 406.0#
            .ka(30) = 395.0#
            .ka(31) = 333.0#
            .ka(32) = 216.0#
            .ka(33) = 112.0#
            .ka(34) = 40.0#
            .ka(35) = 7.5#
        End With
    End Sub ' setkO3

**71

	We run it.  It makes 136 cycles before blowing up.  NaN ("Not A Number") everywhere!
	We have run into something called numerical instability.  A problem with REMs and RCMs back to Manabe and Strickler.  Things are changing too fast for the program to handle.  There's a simple fix:  Decrease the time step.
	We've been going a day at a time.  Let's make it a fraction of one day.
	We will need another global 'run' variable:

    Public timeStep As Double = 0.3#    ' In days

	One modification to the code in the big loop in Main:

            days += timeStep

	And one modification to heat:

        Dim factor As Double = timeStep * 86400.0# * area

	Then we run the program and get:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.183        Fv         62.2
--  -----  -----  -----  -----          TOA    340.3        Fi        277.9
 1   9934  24089  -0.25  169.5          F      278.1        sum       340.1
 2   3279  17483   5.23  167.9          Te     264.6     K  error     -0.051  %
 3   2133  14776   6.63  182.0
 4   1628  12896   7.75  194.5          cycle  1574         res     0.00010   K
 5   1341  11411   8.74  206.0          days   472.2     d  time     3' 32"
 6   1155  10163   9.58  216.9          H2O    0.00e+000 kg
 7   1025   9072  10.28  227.4
 8    927   8096  10.82  237.4          Tstrat 178.5     K  TOA ht    29056   m
 9    851   7207  11.23  247.0          Ttrop  274.0     K  Tpause    12082   m
10    790   6386  11.54  256.2          Tsurf  337.0     K  ground        0   m
11    740   5621  11.77  265.1
12    697   4903  11.94  273.5          r2
13    660   4225  12.05  281.6          rms              K
14    628   3581  12.13  289.4
15    600   2967  12.17  296.8
16    575   2380  12.20  304.0
17    552   1816  12.22  310.8
18    531   1275  12.28  317.5
19    513    753  12.57  323.9
20    496    248  27.43  330.2
21      1     -1   0.00  337.0

	We now have a relatively sexy tool:  a non-gray radiative equilibrium model.  Unfortunately, at 0.3 simulated days (7.2 simulated hours) per time step, it takes a long time to run--472.2 simulated days and 3 minutes 32 seconds in real time.  Interestingly, Manabe and Strickler used an 8-hour time step back in 1964.

**72

	Gradually, we add improvements along the way.  The water vapor burden is trivial.  Just add a bit of code to display:

        H2Okg = 0.0#
        For L As Integer = 1 To N
            H2Okg += layer(L).SM(H2O)
        Next
        H2Okg *= area

	Which reveals that we have...  3.41 x 1015 kg of water vapor in our model atmosphere.  Wait, what?  A hot climate like this, and we've only got a quarter of the 1.27 x 1016 kg we expect?
	Uh oh, a logic error.  We're only calling distrib and setRadAir once, before the main loop starts.  But as temperature changes, so will the saturation vapor pressure and the water vapor content of the air, and thus the radiative properties of each layer.  So we need to call distrib and setRadAir, in that order, in the main loop in the main program, just before the call to transfer.

**73

	And we run it:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.175        Fv         59.6
--  -----  -----  -----  -----          TOA    340.3        Fi        280.6
 1  12040  23520  -0.90  205.5          F      280.7        sum       340.2
 2   3874  15563   0.27  198.3          Te     265.2     K  error     -0.030  %
 3   2334  12459  -0.41  199.2
 4   1660  10461  -0.71  198.4          cycle  1748         res     0.00010   K
 5   1285   8989  -0.53  197.3          days   524.4     d  time     3' 51"
 6   1048   7823  -0.11  196.7          H2O    7.99e+014 kg
 7    886   6856   0.39  196.6
 8    769   6028   0.88  196.9          Tstrat 200.3     K  TOA ht    29540   m
 9    681   5303   1.35  197.6          Ttrop  210.1     K  Tpause     9631   m
10    612   4657   1.84  198.4          Tsurf  306.3     K  ground        0   m
11    557   4072   2.38  199.5
12    511   3538   3.06  200.8          r2
13    474   3045   4.00  202.3          rms              K
14    443   2587   5.45  204.1
15    417   2157   8.03  206.5
16    396   1750  13.36  209.7
17    382   1361  25.74  214.9
18    376    982  49.23  224.7
19    385    602  72.16  243.4
20    409    204 166.89  272.1
21      1     -1   0.00  306.3

	Yick.  Even less water vapor, terribly cold atmosphere except right at the ground.  What's going on here?
	We've added water-vapor feedback to our model.  Under present conditions, cold air and little water vapor feed back on each other, and we wind up quite chilly.
	Perhaps this is because we're still not moving heat from the ground to the air by convection.  Maybe it's time to add that process.

**74

	Back to the lapse rate.  In an atmosphere of a given composition, the "dry" adiabatic lapse rate is simply

	Γ = -g / cp										|39

where g is local gravity and cp the specific heat capacity at constant pressure.  We're using g = 9.80665 m s-2 and cp = 1006 J K-1 kg-1, so Γ = 
-0.009748 K m-1 or about 9.7 K km-1.
	But this is in a dry atmosphere.  In an atmosphere with a big volatile content--like ours, with its 0.4% water vapor by volume--transfer of latent heat due to phase changes alters Γ greatly.  I won't bother with the very complicated equation for the saturated lapse rate, but suffice to say, in the real atmosphere Γ ranges from 4.75 K km-1 near the ground to the dry rate near the tropopause.  A good average for the whole troposphere is about 6.5 K km-1.
	In the real atmosphere, if lapse rates get higher than the environmental rate, convection sets in.  Hot packets of air rise, cool packets sink, and heat gets mixed.  The environmental lapse rate is restored.  One big effect of this is that a lot of heat moves from the ground to the air--on global annual average, 97 watts per square meter.  Of this, 80 is "latent" heat transfer (phase changes), and 17 is "sensible" heat transfer--pure convection plus conduction.
	Therefore, the simplest possible "convective adjustment" for a column atmosphere model is to enforce the environmental lapse rate.  If any layer winds up with too high a lapse rate (and look at the ridiculous rates in the lower layers of our model so far), mix heats with the next lowest layer until you reach the desired level.
	But after you mix the heat from layers 2 and 3, 3 and 4 might be left with too high a rate again.  So you have to keep cycling through the layers until they all have a low enough lapse rate.  Then you've successfully applied a "convective adjustment."
	Once we add this to our model, it will no longer be a radiative equilibrium model.  It will be an honest to God radiative-convective model.

**75

	Here's the convective adjustment routine:

    ' adjust applies the convective adjustment.
    Sub adjust()
        Const crit = 6.5001#

        Do
            For L As Integer = N To 1 Step -1
                If layer(L).gamma > crit Then
                    mix(L)
                End If
            Next

            lapse()

            For L As Integer = 1 To N
                If layer(L).gamma > crit Then
                    Continue Do
                End If
            Next

            Exit Do
        Loop
    End Sub ' adjust

	Note:  I use 6.5001 for the critical level of gamma instead of just 6.5.  This fixes a problem with how computers store floating-point numbers.  Not all numbers can be stored exactly (think of 1/3).  This is kind of a clumsy fix, but hey, it works.

**76

	Here's the heat-mixing routine.  It looks fairly simple.  It ain't.  It took me a long, long time to work out a good algorithm for this.  I won't bore you with the page after page of calculations behind the derivation.

    ' mix redistributes heat among two adjacent layers.
    Sub mix(ByVal L As Integer)
        Const Factor As Double = -0.0065#

        Dim M As Integer = L + 1
        Dim dz As Double = layer(L).z - layer(M).z
        Dim dT As Double = Factor * dz
        Dim IL As Double = 1.0# / layer(L).cpm
        Dim IM As Double = 1.0# / layer(M).cpm

        Dim dH As Double = (dT - layer(L).temp + layer(M).temp)
        dH /= (IL + IM)

        layer(L).temp += dH * IL
        layer(M).temp -= dH * IM
    End Sub ' mix

	Lastly, we need to call adjust from Main, right after the call to heat.

**77

	Results:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.175        Fv         59.5
--  -----  -----  -----  -----          TOA    340.3        Fi        280.5
 1  12434  26938  -0.39  212.2          F      280.7        sum       340.1
 2   4082  18680   0.70  209.0          Te     265.3     K  error     -0.051  %
 3   2476  15401   3.13  211.3
 4   1825  13250   6.50  218.1          cycle  1220         res     0.00010   K
 5   1490  11592   6.50  228.8          days   366.0     d  time     2' 46"
 6   1267  10214   6.50  237.8          H2O    2.51e+016 kg
 7   1107   9028   6.50  245.5
 8    986   7982   6.50  252.3          Tstrat 212.7     K  TOA ht    33155   m
 9    891   7043   6.50  258.4          Ttrop  272.4     K  Tpause    12337   m
10    814   6191   6.50  263.9          Tsurf  304.2     K  ground        0   m
11    751   5409   6.50  269.0
12    697   4685   6.50  273.7          r2
13    652   4010   6.50  278.1          rms              K
14    612   3378   6.50  282.2
15    578   2783   6.50  286.1
16    548   2220   6.50  289.7
17    521   1686   6.50  293.2
18    496   1178   6.50  296.5
19    475    692   6.50  299.7
20    455    227   6.50  302.7
21      1     -1   0.00  304.2

	Now we've got a hotter atmosphere, twice as much more water vapor than we need, and a cooler surface.  304 K.  And we have a working radiative-convective model of Earth's climate system.
	The only major component we're missing now is clouds.  But before we add clouds--quite a complex undertaking, I'm afraid--let's complete the display.  We have two fields for statistical checks.  But what are we checking against?
	The US Standard Atmosphere of 1976.  We will create a table of what the temperatures "should be" at each pressure level according to the USSA-76.  Then we can compare our results and see how close we've gotten.

**78

	Here are the temperatures in the USSA-76 corresponding to the pressure levels in our model:

L#
	p (Pa)
	TUSSA (K)

 1	2474	217.24	
 2	7423	216.66	
 3	12372	216.66	
 4	17321	216.66	
 5	22270	216.75	
 6	27219	224.38	
 7	32167	231.62	
 8	37116	238.02	
 9	42065	243.76	
10	47014	248.99	
11	51963	253.76	
12	56912	258.20	
13	61860	262.34	
14	66809	266.19	
15	71758	269.83	
16	76707	273.28	
17	81656	276.54	
18	86605	279.66	
19	91553	282.64	
20	96502	285.49	
21
	106351
	288.16

**79

	Consider two arrays of real numbers, X() and Y(), each with N elements.  Then, the following routines give the correlation coefficient and the RMS error, respectively.

    ' correl finds the Pearson product-moment correlation coefficient.
    Function correl(ByVal x() As Double, ByVal y() As Double, ByVal N As Integer)
        Dim sumx As Double = 0.0#
        Dim sumy As Double = 0.0#
        Dim sumx2 As Double = 0.0#
        Dim sumy2 As Double = 0.0#
        Dim sumxy As Double = 0.0#
        Dim xN As Double = CDbl(N)

        For i As Integer = 1 To N
            sumx += x(i)
            sumy += y(i)

            sumx2 += x(i) ^ 2
            sumy2 += y(i) ^ 2
            sumxy += x(i) * y(i)
        Next

        Dim meanx As Double = sumx / xN
        Dim meany As Double = sumy / xN
        Dim ssxx As Double = sumx2 - xN * meanx ^ 2
        Dim ssyy As Double = sumy2 - xN * meany ^ 2
        Dim ssxy As Double = sumxy - xN * meanx * meany

        Return ssxy / Sqrt(ssxx * ssyy)
    End Function ' correl

    ' rmse returns the root-mean-squared error.
    Function rmse(ByVal x() As Double, ByVal y() As Double, ByVal N As Integer) As Double
        Dim sum As Double = 0.0#
        Dim xN As Double = CDbl(N)

        For i As Integer = 1 To N
            sum += (x(i) - y(i)) ^ 2
        Next

        Return Sqrt(sum / xN)
    End Function ' rmse

**80

	This routine then shows the statistical fit to the 1976 US Standard Atmosphere:

    ' showStats shows the fit to the USSA-76.
    Sub showStats()
        Const colB As Integer = 48

        Dim Tm(All) As Double
        Dim Tu() As Double = {0.0#, _
        217.24#, 216.66#, 216.66#, 216.66#, 216.75#, 224.38#, 231.62#, _
        238.02#, 243.76#, 248.99#, 253.76#, 258.2#, 262.34#, 266.19#, _
        269.83#, 273.28#, 276.54#, 279.66#, 282.64#, 285.49#, 288.16#}

        For L As Integer = 1 To ground
            With layer(L)
                Tm(L) = .temp
            End With
        Next

        r2 = correl(Tm, Tu, All)
        rms = rmse(Tm, Tu, All)

        put(14, colB, r2, 6, 4)
        put(15, colB, rms, 6, 3)
    End Sub ' showStats

**81

	After the end of the main loop in the main program, call display one last time, then showStats, then pause.  Run the program, and the final screen reads:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.175        Fv         59.5
--  -----  -----  -----  -----          TOA    340.3        Fi        280.5
 1  12434  26938  -0.39  212.2          F      280.7        sum       340.1
 2   4082  18680   0.70  209.0          Te     265.3     K  error     -0.051  %
 3   2476  15401   3.13  211.3
 4   1825  13250   6.50  218.1          cycle  1220         res     0.00010   K
 5   1490  11592   6.50  228.8          days   366.0     d  time     2' 42"
 6   1267  10214   6.50  237.8          H2O    2.51e+016 kg
 7   1107   9028   6.50  245.5
 8    986   7982   6.50  252.3          Tstrat 212.7     K  TOA ht    33155   m
 9    891   7043   6.50  258.4          Ttrop  272.4     K  Tpause    12337   m
10    814   6191   6.50  263.9          Tsurf  304.2     K  ground        0   m
11    751   5409   6.50  269.0
12    697   4685   6.50  273.7          r2     0.9860
13    652   4010   6.50  278.1          rms    14.134    K
14    612   3378   6.50  282.2
15    578   2783   6.50  286.1
16    548   2220   6.50  289.7
17    521   1686   6.50  293.2
18    496   1178   6.50  296.5
19    475    692   6.50  299.7
20    455    227   6.50  302.7
21      1     -1   0.00  304.2

	The correlation is surprisingly good.  r2 is 99%.  RMS error is not so good, a full 14 K.
	It is not too far off to say that r2 measures how well the shape of each curve matches, and RMS error shows how close their location and size are.  Both measures should be good (r2 as close to 1 as possible; RMSE as close to 0 as possible) to declare that we have a good fit.  The trick now is to improve the fit, never by fudging, but always by improving how realistic the physics are.

**82

	There are almost as many cloud schemes in the literature as band schemes.  Manabe's team (Manabe and Strickler 1964, Manabe and Wetherald 1967) used the following scheme:

Clouds	
Av	
Rv	
z	
n	
Ai

High	
0.005	
0.205	
10.0	
0.228	
0.5	
Mid	0.020	0.48	 4.1	0.090	1.0	
Low	0.035	0.69	 2.2	0.313	1.0

	Here Av is shortwave absorptance, Rv shortwave reflectivity, z mean altitude in km, n "cloud amount" (sky fraction covered), and Ai longwave absorptivity.  Meanwhile, Kiehl and Trenberth (1997) used this scheme:

cloud	
km	
n	
Ai

high	
10.5	
0.20	
0.6	
mid	 5.5	0.06	1.0	
low	 1.5	0.49	1.0

	Rossow and Schiffer (1999) had the advantage of accumulated data from the International Satellite Cloud Climatology Project.  They estimate high clouds cover 21.6%, middle 19.3%, low 26.5% of the sky, respectively.
	Note that, assuming random overlap, the Manabe scheme has a net global cloud cover of 52%, Kiehl and Trenberth 62%, and Rossow and Schiffer 54%.
	I'll go with what seems the best data to me--cloud coverage from Rossow and Schiffer, longwave scheme from Kiehl and Trenberth, shortwave scheme (because it's the only one available) from Manabe et al.:

model	
Av	
Rv	
Tv	
Ai	
Ri	
Ti

High	
0.005	
0.205	
0.790	
0.6	
0.0	
0.4	
Mid	0.020	0.480	0.500	1.0	0.0	0.0	
Low	0.035	0.690	0.275	1.0	0.0	0.0

**83

	For altitude, I'll just take the averages for Manabe et al., and Kiehl and Trenberth:

[z]	
MW67	
KT97	
Mean

High	
10.0	
10.5	
10.25	
Mid	 4.1	 5.5	 4.80	
Low	 2.2	 1.5	 1.85

	The pressures for the cloud altitudes, and the closest match for those figures in our layer scheme, are as follows:

z (m)	
p (Pa), USSA-76	
cloud	
layer

10250	
25520	
high	
 6	
 4800	55530	middle	12	
 1850	81019	low	17

**84

	As with gases, we can find A, R and T in each band for cloud particles (Acld, Rcld, and Tcld, actually).  But it is a much more complicated process, as scattering in the shortwave range is significant and we have to use Mie scattering theory.  So for the moment, I'll just go with the figures as given above.  We can calculate our own cloud parameters later.

**85

	To accommodate clouds, we need to make some changes.  In setAbsorbers, add two lines of code:

        set1absorber(ice, 2030.0#, 18.0153#, "ice", 0.0#, 0.0#)
        set1absorber(wtr, 1860.0#, 18.0153#, "wtr", 0.0#, 0.0#)

	Here's a cloud setup routine:

    ' setClouds sets up cloud layers
    Sub setClouds()
        For L As Integer = 1 To N
            With layer(L)
                Select Case L
                    Case cldHi : .fcld = 0.216#
                    Case cldMid : .fcld = 0.193#
                    Case cldLo : .fcld = 0.265#
                    Case Else : .fcld = 0.0#
                End Select

                .fclr = 1.0# - .fcld

                For b As Integer = swlo To swhi
                    Select Case L
                        Case cldHi
                            .Acld(b) = 0.005#
                            .Rcld(b) = 0.205#
                            .Tcld(b) = 0.79#
                        Case cldMid
                            .Acld(b) = 0.02#
                            .Rcld(b) = 0.48#
                            .Tcld(b) = 0.5#
                        Case cldLo
                            .Acld(b) = 0.035#
                            .Rcld(b) = 0.69#
                            .Tcld(b) = 0.275#
                        Case Else
                            .Acld(b) = 0.0#
                            .Rcld(b) = 0.0#
                            .Tcld(b) = 1.0#
                    End Select
                Next

                For b As Integer = lwlo To lwhi
                    Select Case L
                        Case cldHi
                            .Acld(b) = 0.6#
                            .Rcld(b) = 0.0#
                            .Tcld(b) = 0.4#
                        Case cldMid, cldLo
                            .Acld(b) = 1.0#
                            .Rcld(b) = 0.0#
                            .Tcld(b) = 0.0#
                        Case Else
                            .Acld(b) = 0.0#
                            .Rcld(b) = 0.0#
                            .Tcld(b) = 1.0#
                    End Select
                Next
            End With
        Next
    End Sub ' setClouds

	Call this from setup, right after the call to distrib.

**86

	Last, we rewrite setRadAir to take account of our clouds:

    ' setRadAir sets air layer radiative properties.
    Sub setRadAir()
        For L = 1 To N
            With layer(L)
                For b As Integer = 1 To Nbands
                    Dim Tch4 As Double = Exp(-tau(CH4, L, b))
                    Dim Tco2 As Double = Exp(-tau(CO2, L, b))
                    Dim Th2o As Double = Exp(-tau(H2O, L, b))
                    Dim To3 As Double = Exp(-tau(O3, L, b))
                    Dim Tclear As Double = Tch4 * Tco2 * Th2o * To3

                    Select Case L
                        Case cldHi, cldMid, cldLo
                            .T(b) = .fclr * Tclear + .fcld * .Tcld(b) * Tclear
                        Case Else
                            .T(b) = Tclear
                    End Select

                    Dim Rclear As Double = 0.0#
                    If b <= swhi Then
                        Rclear = 0.06# * .mFrac
                    Else
                        Rclear = 0.0#
                    End If

                    Select Case L
                        Case cldHi, cldMid, cldLo
                            .R(b) = .fclr * Rclear + .fcld * .Rcld(b)
                        Case Else
                            .R(b) = Rclear
                    End Select

                    If (.T(b) + .R(b)) >= 1.0 Then
                        .A(b) = 0.0#
                        .T(b) = 1.0# - .R(b)
                    Else
                        .A(b) = 1.0# - .R(b) - .T(b)
                    End If
                Next b
            End With
        Next
    End Sub ' setRadAir

**87

	Run that, and we get:

L#     Av     Rv     Tv       Ai     Ri     Ti
--    -----  -----  -----    -----  -----  -----
 1    0.497  0.003  0.500    0.170  0.000  0.830
 2    0.505  0.003  0.492    0.304  0.000  0.696
 3    0.470  0.003  0.527    0.321  0.000  0.679
 4    0.451  0.003  0.546    0.334  0.000  0.666
 5    0.431  0.003  0.566    0.343  0.000  0.657
 6    0.387  0.047  0.567    0.433  0.000  0.567
 7    0.386  0.003  0.611    0.355  0.000  0.645
 8    0.375  0.003  0.622    0.361  0.000  0.639
 9    0.371  0.003  0.626    0.366  0.000  0.634
10    0.367  0.003  0.630    0.370  0.000  0.630
11    0.367  0.003  0.630    0.375  0.000  0.625
12    0.334  0.095  0.571    0.498  0.000  0.502
13    0.366  0.003  0.631    0.382  0.000  0.618
14    0.367  0.003  0.630    0.385  0.000  0.615
15    0.368  0.003  0.629    0.388  0.000  0.612
16    0.369  0.003  0.628    0.391  0.000  0.609
17    0.307  0.185  0.508    0.554  0.000  0.446
18    0.368  0.003  0.629    0.396  0.000  0.604
19    0.368  0.003  0.629    0.398  0.000  0.602
20    0.367  0.003  0.630    0.400  0.000  0.600
21    0.860  0.140  0.000    0.996  0.004  0.000

	Followed, eventually, by:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.337        Fv        114.6
--  -----  -----  -----  -----          TOA    340.3        Fi        225.4
 1  13887  27235  -1.97  237.0          F      225.6        sum       340.1
 2   4280  18151  -0.90  219.1          Te     251.2     K  error     -0.05   %
 3   2532  14746  -0.34  216.1
 4   1802  12579   1.30  215.3          cycle  1305         res     0.00010   K
 5   1415  10970   5.78  217.4          days   391.5     d  time     2' 57"
 6   1198   9663   6.50  225.0          H2O    8.84e+015 kg
 7   1047   8540   6.50  232.3
 8    932   7551   6.50  238.7          Tstrat 221.9     K  TOA ht    34178   m
 9    843   6663   6.50  244.5          Ttrop  257.7     K  Tpause    11677   m
10    770   5857   6.50  249.7          Tsurf  287.8     K  ground        0   m
11    710   5117   6.50  254.5
12    660   4432   6.50  259.0          r2     0.9867
13    617   3794   6.50  263.1          rms    4.411     K
14    579   3196   6.50  267.0
15    547   2633   6.50  270.7
16    518   2100   6.50  274.1
17    493   1595   6.50  277.4
18    470   1114   6.50  280.5
19    449    655   6.50  283.5
20    430    215   6.50  286.4
21      1     -1   0.00  287.8

	Wow.  r2 is still about 99%, but rms error has dropped to 4.4 K.  The albedo is much closer to normal, though still a bit outside the accurate range of 0.28-0.33 from satellite observations and modeling.  Adding clouds helps a lot.
	Now to do it from physics.

**88

	Let's discuss scattering.  For very small particles, like atoms and molecules, the theory we want would be "Rayleigh scattering," which is pretty much even in all directions.
	For clouds particles, which are larger, we need "Mie scattering" theory--and here, direction is usually more forward than backward.  The contrast between Rayleigh and Mie scattering is why clear skies are blue and cloudy or polluted skies are white.
	Remember those extra terms we tacked onto the absorber structure before:  g, the "asymmetry factor," and ω, the "single-scattering albedo:"

	    ks
	ω = ——										|40
	    ke

	To find out the radiative properties of a cloud, we need to define yet more terms:

	Γ = 2 (1 - ω)0.5 (1 - ω g)0.5							|41

	     (1 - ω g)0.5 - (1 - ω)0.5
	R∞ = ————————————————————————						|42
	     (1 - ω g)0.5 + (1 - ω)0.5

	Γ is hard to explain other than by its mathematical definition.  But R∞ is easy--it's the reflectance a cloud of infinite depth would have.
	With these terms, we can now find R and T for a cloud with no absorption (R and T only):

	      (1 - g) τe
	R = ——————————————								|43
	    1 + (1 - g) τe

	         1
	T = ——————————————								|44
	    1 + (1 - g) τe

	For clouds with some absorption (A, R, and T, the more common case):

	     R∞ [exp(Γ τe) - exp(-Γ τe)]
	 R = ———————————————————————————						|45
	      exp(Γ τe) - R∞2 exp(-Γ τe)

	               1 - R∞2
	 T = ———————————————————————————						|46
	      exp(Γ τe) - R∞2 exp(-Γ τe)

	In each case, the denominator--the bottom half of the fraction--is always the same.  For this discussion I relied heavily on the work of Petty (2006).

**89

	To get any of this stuff, we first need g and ω.  The coefficients of absorption, scattering, and extinction are related by simple algebra:

	ke = ka + ks									|47

	Single-scattering albedo is also pretty simple:

	    ks     ks
	ω = —— = ———————									|48
	    ke   ka + ks

	g is a little harder.  The asymmetry factor ranges from -1 (total backscattering) to 1 (total forward-scattering), and to find exact values for a given particle at a given wavelength, you need a "phase function" to relate everything.  Phase functions can be pretty complicated.
	But we can step around it a bit.  Most of these particle properties can be "parameterized" based on particle size--specifically, the "effective radius" re.  Cloud particle sizes are typically measured in microns.  Chou et al. (1998) assumed 60-75 microns for ice particles and 8-16 for liquid water droplets.  Kiehl and Trenberth (1997) used 10 microns for everything.  The UKMO met office (and its HadGCM climate model) uses 30 microns for ice, 7 for water.  Many field workers have gone out and actually collected particles from balloons, airplanes, helicopters, etc. and measured them, cataloging typical particle sizes for all kinds of clouds.
	A real cloud, of course, almost never has a uniform particle size--there are big drops and little drops, depending on local conditions.  Most analyses assume a "Gaussian" particle distribution.  That's why re is the "effective" particle radius, a mean value that can be substituted for the actual distribution without too much inaccuracy.

**90

	There are lots of re based parameterizations for particle radiative properties.  Because of its combination of relative simplicity and high accuracy, one of the most commonly used is that of Chou et al. 1998:

	ka = a0 + a1 / re									|49
	(1 - ω) = b0 + b1 re + b2 re2							|50
	g = c0 + c1 re + c2 re2								|51

	Here a0 and a1 are constants, as are b0, b1, and b2, and also c0, c1, and c2.  re is the effective radius of the particle.  Given ke and ω, it's simple (see the last section) to find ka and ks.
	The constants are different in different wavelength ranges.  Here are the bands used for this scheme:

	Band I:	0.18 -  0.70 microns
	Band II:	0.70 -  1.22
	Band III:	1.22 -  2.27
	Band IV:	2.27 - 10.00

	To use these, we'll have to incorporate these band boundaries into our band scheme.  The Chou et al. model is for sunlight; note that it doesn't extend past 10 μ.  But we know ahead of time that there's almost no scattering in the infrared under Earth-atmosphere conditions, so all we really need to know is ka.  For that, we can use the UKMO figures:  in the thermal infrared, ka = 65 m2 kg-1 for ice clouds, 130 m2 kg-1 for water clouds.
	Note:  ka in Chou et al. 1998 is in the peculiar units of square meters per gram.  We need to multiply by 1,000 to get SI units.

**91

	To use the Chou et al. cloud parameterizations, we will need no fewer than five new subroutines, plus changes to setClouds.  But first, we need a new data structure:

    Public Structure ChouBandType
        Dim lo As Double
        Dim hi As Double

        Dim g As Double
        Dim ka As Double
        Dim ke As Double
        Dim ks As Double
        Dim omega As Double
    End Structure ' ChouBandType

	And two new global variables:

    Public iceBand(NChou) As ChouBandType
    Public wtrBand(NChou) As ChouBandType

**92

	These first two routines find the optical properties of cloud particles.  The coefficients are from the original article.

    ' setChouIce finds ice particle parameters.
    Sub setChouIce()
        Dim a0() As Double = {0.0#, 0.000333#, 0.000333#, 0.000333#, 0.000333#}
        Dim a1() As Double = {0.0#, 2.52#, 2.52#, 2.52#, 2.52#}

        Dim b0() As Double = {0.0#, 0.0#, -0.0000026#, 0.00215#, 0.0894#}
        Dim b1() As Double = {0.0#, 0.0#, 0.00000746#, 0.000737#, 0.00299#}
        Dim b2() As Double = {0.0#, 0.0#, 0.0#, -0.00000134#, -0.0000104#}

        Dim c0() As Double = {0.0#, 0.746#, 0.749#, 0.761#, 0.841#}
        Dim c1() As Double = {0.0#, 0.000105#, 0.0012#, 0.00142#, 0.00126#}
        Dim c2() As Double = {0.0#, -0.00000264#, -0.00000367#, -0.00000396#, _
                             -0.00000385#}

        Dim re As Double = ri
        Dim re2 As Double = re ^ 2

        For cb As Integer = 1 To 4
            With iceBand(cb)
                .ke = 1000.0# * (a0(cb) + a1(cb) / re)
                Dim comp As Double = b0(cb) + b1(cb) * re + b2(cb) * re2
                .omega = 1.0# - comp
                .ks = .omega * .ke
                .ka = .ke - .ks
                .g = c0(cb) + c1(cb) * re + c2(cb) * re2
            End With
        Next

        With iceBand(5)
            .ke = 65.0#
            .omega = 0.0#
            .ks = 0.0#
            .ka = 65.0#
            .g = 0.92#
        End With
    End Sub ' setChouIice

    ' setChouWtr finds water droplet parameters.
    Sub setChouWtr()
        Dim a0() As Double = {0.0#, -0.00659#, -0.0101#, -0.0166#, -0.0339#}
        Dim a1() As Double = {0.0#, 1.65#, 1.72#, 1.85#, 2.16#}

        Dim b0() As Double = {0.0#, 0.0#, 0.0000000715#, 0.000888#, 0.0179#}
        Dim b1() As Double = {0.0#, 0.0#, 0.00000845#, 0.000888#, 0.0179#}
        Dim b2() As Double = {0.0#, 0.0#, -0.00000000415#, -0.0000065#, _
                             -0.000369#}

        Dim c0() As Double = {0.0#, 0.826#, 0.749#, 0.745#, 0.835#}
        Dim c1() As Double = {0.0#, 0.00529#, 0.00832#, 0.0137#, 0.00257#}
        Dim c2() As Double = {0.0#, -0.00000149#, -0.000233#, -0.000382#, _
                             -0.0000552#}

        Dim re As Double = rw
        Dim re2 As Double = re ^ 2

        For cb As Integer = 1 To 4
            With wtrBand(cb)
                .ke = 1000.0# * (a0(cb) + a1(cb) / re)
                Dim comp As Double = b0(cb) + b1(cb) * re + b2(cb) * re2
                .omega = 1.0# - comp
                .ks = .omega * .ke
                .ka = .ke - .ks
                .g = c0(cb) + c1(cb) * re + c2(cb) * re2
            End With
        Next

        With wtrBand(5)
            .ke = 130.0#
            .omega = 0.0#
            .ks = 0.0#
            .ka = 130.0#
            .g = 0.85#
        End With
    End Sub ' setChouWtr

**93

	Next comes a routine to find R and T for clouds, now that we know the particle properties:

    ' setCloudRT finds cloud radiative properties.
    Sub setCloudRT(ByVal medium As Integer, ByVal L As Integer, ByVal b As Integer)
        Dim denom As Double = 0.0#

        With absorber(medium)
            Dim taue As Double = .ke(b) * layer(L).SM(medium)

            If .om(b) = 1.0# Then
                Dim gt As Double = (1.0# - .g(b)) * taue

                denom = 1.0# + gt
                layer(L).Rcld(b) = gt / denom
                layer(L).Tcld(b) = 1.0# / denom
            Else
                Dim o As Double = Sqrt(1.0# - .om(b))
                Dim og As Double = Sqrt(1.0# - .om(b) * .g(b))
                Dim Gamma As Double = 2.0# * o * og
                Dim Rinf As Double = (og - o) / (og + o)
                Dim e1 As Double = Exp(Gamma * taue)
                Dim e2 As Double = Exp(-Gamma * taue)

                denom = e1 - Rinf ^ 2 * e2
                layer(L).Rcld(b) = Rinf * (e1 - e2) / denom
                layer(L).Tcld(b) = (1.0# - Rinf ^ 2) / denom
            End If
        End With
    End Sub ' setCloudRT

	Note that we use .om for the single-scattering albedo in absorbers, rather than the .omega we use in Chou bands.  This has no particular significance other than inconsistency on the part of the author (me).

**94

	Just as we need "setk" routines for the greenhouse gases, we need them for cloud particles.  We've found properties for particles in Chou bands, now we need them for absorbers in our own band scheme.  Here's the code:

    ' setkIce sets ice radiative properties for each band.
    Sub setkIce()
        With absorber(ice)
            For b As Integer = 1 To Nbands
                Dim cIndex As Integer = 0
                Dim midPoint As Double = 0.5# * (band(b).lo + band(b).hi)

                If midPoint < 0.7# Then
                    cIndex = 1
                ElseIf midPoint >= 0.7# And midPoint < 1.22# Then
                    cIndex = 2
                ElseIf midPoint >= 1.22 And midPoint < 2.27# Then
                    cIndex = 3
                ElseIf midPoint >= 2.27# And midPoint < 4.0# Then
                    cIndex = 4
                Else
                    cIndex = 5
                End If

                .g(b) = iceBand(cIndex).g
                .ka(b) = iceBand(cIndex).ka
                .ke(b) = iceBand(cIndex).ke
                .ks(b) = iceBand(cIndex).ks
                .om(b) = iceBand(cIndex).omega
            Next
        End With
    End Sub ' setkIce

    ' setkWtr sets water radiative properties for each band.
    Sub setkWtr()
        With absorber(wtr)
            For b As Integer = 1 To Nbands
                Dim cIndex As Integer = 0
                Dim midPoint As Double = 0.5# * (band(b).lo + band(b).hi)

                If midPoint < 0.7# Then
                    cIndex = 1
                ElseIf midPoint >= 0.7# And midPoint < 1.22# Then
                    cIndex = 2
                ElseIf midPoint >= 1.22 And midPoint < 2.27# Then
                    cIndex = 3
                ElseIf midPoint >= 2.27# And midPoint < 4.0# Then
                    cIndex = 4
                Else
                    cIndex = 5
                End If

                .g(b) = wtrBand(cIndex).g
                .ka(b) = wtrBand(cIndex).ka
                .ke(b) = wtrBand(cIndex).ke
                .ks(b) = wtrBand(cIndex).ks
                .om(b) = wtrBand(cIndex).omega
            Next
        End With
    End Sub ' setkWtr

	Don't forget to call these from the end of setBands!

**95

	Finally, we need to rewrite setClouds:

    ' setClouds sets up cloud layers
    Sub setClouds()
        setChouIce()
        setChouWtr()
        setkIce()
        setkWtr()

        For L As Integer = 1 To N
            With layer(L)
                Select Case L
                    Case cldHi
                        .fcld = 0.216#
                        .SM(ice) = 0.009#
                        .SM(wtr) = 0.0#
                    Case cldMid
                        .fcld = 0.193#
                        .SM(ice) = 0.0#
                        .SM(wtr) = 0.02#
                    Case cldLo
                        .fcld = 0.265#
                        .SM(ice) = 0.0#
                        .SM(wtr) = 0.036#
                    Case Else
                        .fcld = 0.0#
                        .SM(ice) = 0.0#
                        .SM(wtr) = 0.0#
                End Select

                .fclr = 1.0# - .fcld

                For b As Integer = 1 To Nbands
                    Select Case L
                        Case cldHi
                            setCloudRT(ice, L, b)
                        Case cldMid, cldLo
                            setCloudRT(wtr, L, b)
                        Case Else
                            .Rcld(b) = 0.0#
                            .Tcld(b) = 0.0#
                    End Select

                    .Acld(b) = 1.0# - .Rcld(b) - .Tcld(b)
                    If .Acld(b) < 0.0# Then
                        .Acld(b) = 0.0#
                        .Tcld(b) = 1.0# - .Rcld(b)
                    End If
                Next
            End With
        Next
    End Sub ' setClouds

**96

	We've added a lot of code to do cloud physics.  A lot of code.  All kinds of scope for logic errors.  Before we go any further, it's time to write out more diagnostics.  Here's one for the Chou bands:

    ' showChou shows Chou et al. 1998 band data.
    Sub showChou()
        Clear()
        WriteLine("ice re " & ri.ToString)
        WriteLine("c    g      ka       ks       ke     omega")
        WriteLine("-  -----  -------  -------  -------  -----")

        For cb As Integer = 1 To 5
            Write(cb)

            With iceBand(cb)
                Write(.g.ToString("f3").PadLeft(7))
                Write(.ka.ToString("f3").PadLeft(9))
                Write(.ke.ToString("f3").PadLeft(9))
                Write(.ks.ToString("f3").PadLeft(9))
                Write(.omega.ToString("f3").PadLeft(7))
            End With

            WriteLine()
        Next

        WriteLine()
        WriteLine("wtr re " & rw.ToString)
        WriteLine("c    g      ka       ks       ke     omega")
        WriteLine("-  -----  -------  -------  -------  -----")

        For cb As Integer = 1 To 5
            Write(cb)

            With wtrBand(cb)
                Write(.g.ToString("f3").PadLeft(7))
                Write(.ka.ToString("f3").PadLeft(9))
                Write(.ke.ToString("f3").PadLeft(9))
                Write(.ks.ToString("f3").PadLeft(9))
                Write(.omega.ToString("f3").PadLeft(7))
            End With

            WriteLine()
        Next

        pause()
    End Sub ' showChou

**97

	We call this from diag right after the call to showPhys.  It shows us:

ice re 30
c    g      ka       ke       ks     omega
-  -----  -------  -------  -------  -----
1  0.747    0.000   84.333   84.333  1.000
2  0.782    0.019   84.333   84.314  1.000
3  0.800    1.944   84.333   82.389  0.977
4  0.875   14.315   84.333   70.018  0.830
5  0.920   65.000   65.000    0.000  0.000

wtr re 5
c    g      ka       ke       ks     omega
-  -----  -------  -------  -------  -----
1  0.852    0.000  323.410  323.410  1.000
2  0.785    0.014  333.900  333.886  1.000
3  0.804    1.825  353.400  351.575  0.995
4  0.846   39.083  398.100  359.017  0.902
5  0.850  130.000  130.000    0.000  0.000

	So far, so good.

**97

	Now let's add a diagnostic routine to show what we're getting for the clouds.  We need a new "semigray" routine to find these from the band-by-band values:

    ' semigrayCloud finds two-tone cloud radiative properties.
    Sub semigrayCloud()
        Dim idx() As Integer = {0, cldHi, cldMid, cldLo}

        For i As Integer = 1 To 3
            Dim L As Integer = idx(i)

            With layer(L)
                .Avcld = 0.0#
                .Rvcld = 0.0#
                .Tvcld = 0.0#

                .Aicld = 0.0#
                .Ricld = 0.0#
                .Ticld = 0.0#

                For b As Integer = swlo To swhi
                    .Avcld += .Acld(b)
                    .Rvcld += .Rcld(b)
                    .Tvcld += .Tcld(b)
                Next

                For b As Integer = lwlo To lwhi
                    .Aicld += .Acld(b)
                    .Ricld += .Rcld(b)
                    .Ticld += .Tcld(b)
                Next

                .Avcld /= Nsw
                .Rvcld /= Nsw
                .Tvcld /= Nsw

                .Aicld /= Nlw
                .Ricld /= Nlw
                .Ticld /= Nlw
            End With
        Next
    End Sub ' semigrayCloud

**98

	Finally, here's a cloud diagnostics screen:

    ' showClouds shows cloud radiative properties.
    Sub showClouds()
        semigrayCloud()
        Clear()
        WriteLine("L#    Avcld  Rvcld  Tvcld    Aicld  Ricld  Ticld")
        WriteLine("--    -----  -----  -----    -----  -----  -----")

        Dim idx() As Integer = {0, cldHi, cldMid, cldLo}

        For i As Integer = 1 To 3
            Dim L As Integer = idx(i)
            Write(L.ToString.PadLeft(2))

            With layer(L)
                Write(.Avcld.ToString("f3").PadLeft(9))
                Write(.Rvcld.ToString("f3").PadLeft(7))
                Write(.Tvcld.ToString("f3").PadLeft(7))

                Write(.Aicld.ToString("f3").PadLeft(9))
                Write(.Ricld.ToString("f3").PadLeft(7))
                Write(.Ticld.ToString("f3").PadLeft(7))
            End With

            WriteLine()
        Next

        pause()
    End Sub ' showClouds

	Call this from diag right after the call to showChou.  Run the program, and this routine shows us:

L#    Avcld  Rvcld  Tvcld    Aicld  Ricld  Ticld
--    -----  -----  -----    -----  -----  -----
 6    0.081  0.119  0.801    0.690  0.000  0.310
12    0.240  0.417  0.343    0.994  0.000  0.006
17    0.274  0.505  0.220    1.000  0.000  0.000

**99

	Finally, let's modify the showSM routine to show us cloud mass paths (and while we're at it, better put in methane, since we're now using that as a greenhouse gas):

    ' showSM shows mass paths.
    Sub showSM()
        Clear()
        WriteLine("L#   CH4      CO2      H2O      O3       ice     water ")
        WriteLine("--  -------  -------   -------   -------   -------   -------")

        For L As Integer = 1 To N
            Write(L.ToString.PadLeft(2))

            With layer(L)
                Write(.SM(CH4).ToString("f5").PadLeft(9))
                Write(.SM(CO2).ToString("f3").PadLeft(9))
                Write(.SM(H2O).ToString("f3").PadLeft(9))
                Write(.SM(O3).ToString("f5").PadLeft(9))

                Select Case L
                    Case cldHi
                        Write(.SM(ice).ToString("f3").PadLeft(9))
                    Case cldMid, cldLo
                        Write(.SM(wtr).ToString("f3").PadLeft(18))
                End Select
            End With

            WriteLine()
        Next

        pause()
    End Sub ' showSM

**100

	Running this routine shows us:

L#    CH4      CO2      H2O      O3       ice     water
--  -------  -------  -------  -------  -------  -------
 1  0.00045    0.247    0.067  0.00368
 2  0.00045    0.247    0.271  0.00138
 3  0.00045    0.247    0.312  0.00052
 4  0.00045    0.247    0.329  0.00030
 5  0.00045    0.247    0.339  0.00018
 6  0.00045    0.247    0.345  0.00010    0.009
 7  0.00045    0.247    0.349  0.00007
 8  0.00045    0.247    0.352  0.00005
 9  0.00045    0.247    0.355  0.00004
10  0.00045    0.247    0.356  0.00003
11  0.00045    0.247    0.358  0.00003
12  0.00045    0.247    0.359  0.00003             0.020
13  0.00045    0.247    0.360  0.00003
14  0.00045    0.247    0.361  0.00003
15  0.00045    0.247    0.362  0.00003
16  0.00045    0.247    0.363  0.00003
17  0.00045    0.247    0.363  0.00003             0.036
18  0.00045    0.247    0.364  0.00003
19  0.00045    0.247    0.364  0.00002
20  0.00045    0.247    0.365  0.00002

**101

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.321        Fv        109.2
--  -----  -----  -----  -----          TOA    340.3        Fi        230.9
 1  13903  27406  -1.95  237.3          F      231.1        sum       340.1
 2   4288  18311  -0.85  219.5          Te     252.7     K  error     -0.05   %
 3   2539  14897  -0.18  216.6
 4   1810  12723   1.82  216.3          cycle  1384         res     0.00010   K
 5   1427  11104   6.50  219.2          days   415.2     d  time     3' 10"
 6   1213   9784   6.50  227.8          H2O    1.12e+016 kg
 7   1060   8648   6.50  235.2
 8    944   7646   6.50  241.7          Tstrat 222.4     K  TOA ht    34357   m
 9    853   6747   6.50  247.5          Ttrop  260.9     K  Tpause    11818   m
10    780   5931   6.50  252.8          Tsurf  291.4     K  ground        0   m
11    719   5181   6.50  257.7
12    668   4488   6.50  262.2          r2     0.9890
13    624   3842   6.50  266.4          rms     5.655    K
14    587   3236   6.50  270.3
15    554   2666   6.50  274.0
16    525   2127   6.50  277.6
17    499   1615   6.50  280.9
18    476   1128   6.50  284.0
19    455    663   6.50  287.1
20    436    218   6.50  290.0
21      1     -1   0.00  291.4

	Hmm.  Performance is a little worse this time, despite having more realistic clouds.
	If you improve the simulation's physics, and the fit becomes worse, it means something else you're doing is not accurate enough, and your earlier approximations were covering it up.  Let's stop tinkering around with clouds and see what else we can make more realistic.

**102

	Our whole climate system is a little too hot, but the stratosphere seems the most affected.  What determines the temperature of the stratosphere?
	The balance is set by absorption of UV, mostly by ozone, and emission of IR, mostly by carbon dioxide.  Water vapor is a very minor player because the stratosphere is very dry.  So we're either absorbing too much UV, or radiating too little IR.
	But our band scheme for ozone is quite precise, and our carbon dioxide figures aren't bad, either.  The only possible wrong input would be the level of sunlight itself…
	Which, of course, it is.  I'm using the Planck blackbody relation to set the Solar fraction.  This is only approximately correct.  In reality, the Sun is not quite a Planckian radiator.  It has absorption lines from the hydrogen, helium, and "metals" (respectively 71%, 28%, and 1% of the mass) in its outer layers, and some wavelengths are blocked more than others, and some emit more than you'd expect.  The deviations aren't great...  but they're there.  How can we introduce this into our model?  Will a more realistic Solar spectrum help?

**103

	Fortunately, NASA and other countries' space agencies have already done the work for us.  I summarize what I calculated from a table in Allen (1973).  Here, NPF is a multiplier I call the Non-Planckian Factor:

Band	Low	High	Planck D	NPF	

 1	0.175	0.225	0.00393	0.279	
 2	0.225	0.245	0.00413	0.277	
 3	0.245	0.280	0.01263	0.408	
 4	0.280	0.295	0.00781	0.851	
 5	0.295	0.310	0.00930	0.880	
 6	0.310	0.320	0.00702	0.851	
 7	0.320	0.400	0.07674	0.866	
 8	0.400	0.700	0.36648	1.116	
 9	0.700	1.800	0.43335	1.008	
10	1.800	1.900	0.00964	0.941	
11	1.900	2.000	0.00814	1.114	
12	2.000	2.100	0.00692	0.586	
13	2.100	2.200	0.00592	0.685	
14	2.200	2.500	0.01329	0.915	
15	2.500	2.600	0.00333	0.626	
16	2.600	2.900	0.00774	0.807	
17	2.900	3.000	0.00200	1.039	
18	3.000	3.500	0.00721	0.481	
19	3.500	4.100	0.00494	0.762	
20	4.100	4.500	0.00207	0.566	
21
	4.500
	5.000
	0.00179
	0.817

	You will observe that most of the UV is depleted.  After band #21 it settles down to a very close approximation of the Planck curve.

**104

	So let's rewrite set1band to allow for a non-Planckian Sun:

    ' set1band sets properties for a single band.
    Sub set1band(ByVal b As Integer, ByVal ll As Double, ByVal hh As Double)
        Dim NPfactor() As Double = {0.0#, _
        0.279#, 0.277#, 0.408#, 0.851#, 0.88#, 0.851#, 0.866#, _
        1.116#, 1.008#, 0.941#, 1.114#, 0.586#, 0.685#, 0.915#, _
        0.626#, 0.807#, 1.039#, 0.481#, 0.762#, 0.566#, 0.817#}

        With band(b)
            .lo = ll
            .hi = hh
            .solFrac = D(hh, 5774.0#) - D(ll, 5774.0#)

            If b <= 21 Then
                .solFrac *= NPfactor(b)
            End If
        End With
    End Sub ' set1band

**104

	We then get:

L#  ht(m)  z (m)  lapse  T (K)          Albedo 0.332        Fv        112.9
--  -----  -----  -----  -----          TOA    340.3        Fi        228.8
 1  13323  27057  -1.20  227.4          F      227.3        sum       341.7
 2   4236  18277  -0.55  216.9          Te     251.6     K  error      0.43   %
 3   2520  14900   0.15  215.0
 4   1802  12739   2.61  215.3          cycle  1372         res     0.00010   K
 5   1429  11123   6.50  219.6          days   411.6     d  time     3'  7"
 6   1215   9800   6.50  228.2          H2O    1.16e+016 kg
 7   1062   8662   6.50  235.6
 8    946   7658   6.50  242.1          Tstrat 218.7     K  TOA ht    33718   m
 9    855   6758   6.50  247.9          Ttrop  261.3     K  Tpause    11837   m
10    781   5940   6.50  253.2          Tsurf  291.9     K  ground        0   m
11    720   5190   6.50  258.1
12    669   4495   6.50  262.6          r2     0.9965
13    625   3848   6.50  266.8          rms     4.523    K
14    588   3241   6.50  270.8
15    555   2670   6.50  274.5
16    525   2130   6.50  278.0
17    500   1618   6.50  281.3
18    476   1130   6.50  284.5
19    455    664   6.50  287.5
20    436    218   6.50  290.4
21      1     -1   0.00  291.9

	Which is improved from last time, though still not great.  r2 is very high indeed, but we've only brought rms error down to 4.5 K.  And in addition, the energy conservation error is up by a factor of 8, though still well under 1%.
	How can we improve things more?  A better band scheme would help, especially for water vapor.  Gonima (1992) lists 66 longwave absorption bands for water vapor.  We could improve our Rayleigh scattering scheme, add more absorbers, get more up-to-date pressure and temperature exponents for absorbers, etc., etc.  Use 0.9 for the pressure exponent for water vapor (ask me why), and rms error falls to 3.1 K.
	At any rate, the process should be clear now.  You start with something simple that barely works.  You revise it.  Add more physics.  Lather, rinse, repeat.  Deniers accuse climate scientists of always changing the models.  Duh!  Of course we do, whenever we can.  There are good, honest ways to change the models:  improve the physics, the resolution, the algorithm.  There are bad, dishonest ways to change the models:  use fudge factors, ignore important physics.  Somebody who tried the latter would be called out in peer review (and don't think it doesn't happen).  The honest improvements make it into the journals.
	I've presented a single-column model with 1,599 lines of code--not counting blank lines and comment lines, probably more like 1,500.  A model like this would be state-of-the-art for c. 1970.  The General Circulation Models now in use by climate scientists are incredibly more sophisticated, and are put together not by guys like me working on PCs, but by whole teams working on supercomputers.  They not only continually improve their models, they check each other's work--you can find article after article reviewing the models, criticizing them.  There's even a "model intercomparison project" whose main purpose is to look for flaws and bad assumptions.  The work goes on.  That's what science is all about.

**References

Allen CW 1973.  Astrophysical Quantities, 3rd ed.  London, UK:  The Athlone Press.

Chou M-D, Lee K-T 1996.  Parameterizations for the absorption of solar radiation by water vapor and ozone.  J. Atmos. Sci. 53, 1203-1208.

Chou M-D, Suarez, MJ, Ho C-H, Yan, MM-H, and Lee K-T  1998.  Parameterizations for Cloud Overlapping and Shortwave Single-Scattering Properties for Use in General Circulation and Cloud Ensemble Models.  J. Clim. 11, 202-214.

Essenhigh RJ 2001.  Robert Essenhigh Replies.  Letter, Chemical Innovation 31, 62-64.

Gonima L 1992.  Validation of a radiation model for estimation of longwave net radiation at the surface.  Clim. Res. 2, 55-63.

Hanel RA, Bandeen WL, Conrath BJ 1963.  The Infrared Horizon of the Planet Earth.  NASA Technical Note D-1850.  Greenbelt, MD:  Goddard Space Flight Center.

Hart MH 1978.  The evolution of the atmosphere of the Earth.  Icarus 33, 23-39.

Haurwitz G 1948.  Insolation in relation to cloud type.  J. Meteorol. 5, 110-113.

Hulburt EO 1931.  The Temperature of the Lower Atmosphere of the Earth.  Phys. Rev. 38, 1876-1890.

Kiehl JT, Trenberth KE 1997.  Earth's annual global mean energy budget.  Bull. Am. Meteorol. Soc. 78, 197-208.

Manabe S, Strickler RF 1964.  Thermal equilibrium of the atmosphere with a convective adjustment.  J. Atmos. Sci. 21, 361-385.

Manabe S, Wetherald RT 1967.  Thermal equilibrium of the atmosphere with a given distribution of relative humidity.  J. Atmos. Sci. 24, 241-259.

Petty RW 2006 (2002).  A First Course in Atmospheric Radiation (2nd edition).  Tucson, AZ:  Sundog Press.

Rossow WB, Schiffer RA 1999.  Advances in understanding clouds from ISCCP.  Bull. Am. Meteorol. Soc. 80, 2261-2287.

Trenberth KE, Smith L 2005.  The Mass of the Atmosphere:  A Constraint on Global Analyses.  J. Clim. 18, 864-875.